Time Dependent Geometry in Massive Gravity
Yaghoub Heydarzade Prabir Rudra Behnam Pourhassan Mir Faizal Ahmed Farag Ali Farhad Darabi
physics.gen-ph




a]Yaghoub Heydarzade,
b]Prabir Rudra,
c]Behnam Pourhassan,
d]Mir Faizal,
e]Ahmed Farag Ali,
f]Farhad Darabi

[a] Department of Physics, Azarbaijan Shahid Madani University, Tabriz, 53714-161, Iran
[b] Department of Mathematics, Asutosh College, Kolkata-700 026, India
[c] School of Physics, Damghan University, Damghan, 3671641167, Iran
[d] Irving K. Barber School of Arts and Sciences,  University of British Columbia - Okanagan, Kelowna, BC V1V 1V7, Canada
[d] Department of Physics and Astronomy, University of Lethbridge, Lethbridge,   Alberta, T1K 3M4, Canada  
[e] Netherlands Institute for Advanced Study, Korte Spinhuissteeg 3, 1012 CG Amsterdam, Netherlands
[e] Department of Physics, Faculty of Science, Benha University, Benha, 13518, Egypt
[f] Department of Physics, Azarbaijan Shahid Madani University, Tabriz, 53714-161, Iran


heydarzade@azaruniv.edu
prudra.math@gmail.com
b.pourhassan@du.ac.ir
mirfaizalmir@googlemail.com
ahmed.ali@fsc.bu.edu.eg
f.darabi@azaruniv.edu (Corresponding author)


In this paper, we will analyze a time dependent geometry in a massive theory of gravity.
This will be done by analyzing
 Vaidya spacetime in such a massive theory of gravity.
As gravitational collapse is a time dependent system, we will
analyze it using the Vaidya spacetime in massive gravity. The
Vainshtein  and dRGT mechanisms are used to obtain a ghost free
massive gravity, and construct such time dependent solutions.
Singularities formed, their nature and strength will be studied in
detail. We will also study the thermodynamical aspects of such a
geometry by calculating the important thermodynamical quantities
for such a system, and analyzing the thermodynamical behavior of
such quantities.




Time Dependent Geometry  in  Massive Gravity
    [
    December 30, 2023
============================================



§ INTRODUCTION

The  observations from   type I supernovae indicate that our
universe is in a state of accelerated cosmic expansion <cit.>-<cit.>.
This accelerated cosmic expansion can be explained by a cosmological constant term
in the Einstein equation, and the existence of such a cosmological constant is predicted
from  all quantum field theories. However, the cosmological constant paradigm
suffers from two well known problems as the “cosmological constant problem” and
the “coincidence problem”.  These problems have motivated the people to
do research in the dark energy models <cit.>-<cit.> and the modified theories
of gravity <cit.>. The latter theories are constrained by the solar system tests <cit.>-<cit.>,
where the modifications have to occur only at the infrared limit.
It is possible to obtain an infrared modification of the general relativity
by  making the gravitons massive <cit.>, such that the small graviton mass
does not violate the known experimental bounds.
Even though this  has been done by adding a small Fierz-Pauli mass term to the original action of
general relativity  <cit.>-<cit.>,  there are problems with the zero
mass limit of this theory due to the   force
mediated by the scalar graviton. Furthermore, such a modified theory of gravity
violates the experimental  bounds obtained from solar system experiments, and so it cannot be a physical theory
<cit.>-<cit.>.

It  was possible to resolve these problems
by using the   Vainshtein mechanism, which was based on the  inclusion of
non-linear
terms in the field equation <cit.>-<cit.>. Even though the Vainshtein mechanism
produces the general relativity in suitable limits, it contains higher derivative terms.
These higher derivative terms give rise to negative norm   Boulware-Deser ghosts <cit.>.
This problem can also be resolved for a subclass of massive potentials, as it has been observed
that for such a subclass of massive potential the
Boulware-Deser ghosts do not appear <cit.>-<cit.>.
This has been done using dRGT mechanism, which is
a theory with one dynamical and one fixed metric
<cit.>.  It is interesting to note that a mass term in the gravitational action can
also be generated from the spontaneous breaking of Lorentz symmetry at the cosmological scale.
<cit.>-<cit.>.  Thus, the  massive gravity  might be produced by some interesting
theoretical considerations.

As massive gravity produces interesting deformation of the general relativity, it has
been used to study the behavior of various interesting systems.
The  black holes in Gauss-Bonnet massive gravity have been studied <cit.>-<cit.>, and it has been demonstrated
that the inclusion of mass term produces interesting deformation of these black holes.
The thermodynamics of such black holes has been studied in the extended phase space
<cit.>-<cit.>. It has been demonstrated
that the phase transition of black holes depends on the different parameters used in this
 massive gravity <cit.>.
Cosmological solutions with a well defined   initial values have been constructed in massive gravity <cit.>.
 The initial value constraints in massive gravity  have been used to study the
 spherically symmetric deformations of flat space,
 and it has been demonstrated that there is a physical sector of the theory, where the theory is stable <cit.>.


The massive theory of gravity has also been used to analyze the deformation of AdS spacetime,
and its CFT dual using the AdS/CFT correspondence <cit.>-<cit.>.
The holographic entanglement entropy of a field theory dual to the massive gravity has also
been studied <cit.>. It was observed using this holographic entanglement
entropy that  both  first order phase
transition  and  second order phase transition occur in this system.
The holographic complexity has also been calculated in the massive gravity <cit.>.
The  stability of solution in massive gravity have been studied using
holographic conductivity <cit.>. Thus, massive gravity has  been used
to study interesting physical systems using gravity/gravity duality. This is another
motivation for analyzing solutions in massive gravity. As the massive gravity is interesting
modification of general relativity, we will analyze a time dependent solutions in massive gravity.

The  time dependent deformation of AdS solution has been used to
analyze the time dependent field theories
<cit.>-<cit.>, and it has led some interests to study
such solutions in massive gravity. These solutions are obtained as
deformations of the Vaidya spacetime, which is a time dependent
spherically symmetric spacetime <cit.>-<cit.>. In fact,
a time dependent black hole solution <cit.>, and  a time
dependent solutions in AdS/CFT correspondence <cit.>, have
been studied using massive gravity. The Vaidya spacetime has
already been used to investigate the jet quenching <cit.> of
virtual gluons and thermalization of a strongly-coupled plasma
<cit.>, with a non-zero chemical potential via the
gauge/gravity duality. Thus, Vaidya spacetime can be used to model
interesting physical systems. We would like to point out that
Vaidya spacetime  has also been used to analyze  gravitational
collapse <cit.>-<cit.>. In fact, the gravitational
collapse in Vaidya spacetime has been widely studied in different
scenarios <cit.>-<cit.>. So, the study of gravitational
collapse is an important application to Vaidya spacetime. We would
like to point out that either  black holes or naked singularities
form from such a gravitational collapse. The difference between
these two types of singularities is that black holes are covered
by a boundary known as event horizon beyond which no information
is conveyed to an external observer. But naked singularities are
not covered by any such boundaries. Theoretically, the existence
of naked singularity is important because that would mean that it
is possible to observe the gravitational collapse of an object to
infinite density. The formation of naked singularities in general
relativity have been studied using Vaidya spacetime
<cit.>-<cit.>. As this is an interesting physical
system, and massive gravity is an important modification of
general relativity, we will study the formation of naked
singularities in massive gravity. So, in this paper, we will first
study a time dependent geometry in massive gravity using Vaidya
spacetime. Then we will use this solution  to analyze the
formation of naked singularity in massive gravity.  We will also
study the thermodynamics of such a time dependent solution.



§ VAIDYA SPACETIME IN MASSIVE GRAVITY

In this section, we will study a time dependent geometry using
Vaidya spacetime in  massive gravity. The  four dimensional action of massive gravity is given by
 
    ℐ=∫ d^4x √(-g)[ℛ+
    ℳ^2∑_i^4c_i𝒰_i(g,f)+ℒ_m],

where f is a fixed symmetric tensor and is called the reference metric, c_i are constants, ℳ is the massive gravity
parameter and 𝒰_i are symmetric polynomials of the
eigenvalues of the d× d matrix
𝒦^μ_ν=√(g^μαf_αν)
given by

    𝒰_1=[𝒦],
          𝒰_2=[𝒦]^2 -[𝒦^2],
          𝒰_3=[𝒦]^3 -3[𝒦][𝒦^2]  + 2[𝒦^3],
          𝒰_4=[𝒦]^4 -6[𝒦^2][𝒦]^2  + 8[𝒦^3][𝒦]+3[𝒦^2]^2
    -6[𝒦^4].

The square root in 𝒦 means
(√(A))^μ_ν(√(A))^ν_λ=A^μ_λ and 𝒦=𝒦^μ_μ.
Then, the equation of motion of this action will be

    G_μν+ℳ^2χ_μν=T_μν,

where G_μν is the Einstein tensor and χ _μν is

    χ _μν =   -c_1/2( 𝒰_1g_μν-
    𝒦_μν) -c_2/2( 𝒰_2g_μν-2𝒰_1𝒦_μν+2𝒦_μν^2)
          -c_3/2(𝒰_3g_μν-3𝒰_2𝒦_μν +6𝒰_1𝒦_μν^2-6𝒦_μν^3)
          -c_4/2(𝒰_4g_μν-4𝒰_3𝒦_μν+12𝒰_2𝒦_μν^2-24𝒰_1𝒦 _μν^3+24𝒦_μν^4)

Now, we can investigate a Vaidya metric in the context of massive gravity.
We consider the spatial reference metric, in the basis (t, r, θ, ϕ), as follows
 <cit.>

    f_μν=diag(0, 0, c^2 h_ij),

where h_ij is two dimensional Euclidean metric and c is a
positive constant. The Vaidya metric in the advanced time
coordinate system is given by

    ds^2 =f(t, r)dt^2 +2dtdr+r^2 dΩ_2^2,

 where

    f(t, r)= -(1-m(t,r)/r).

We also consider the supporting total energy-momentum tensor of the
field equation (<ref>) in the following form

    T_μν=T_μν^(n)+T_μν^(m),

where T_μν^(n) and T_μν^(m) are the energy-momentum tensor for the Vaidya null radiation and the energy-momentum
tensor of the perfect fluid
supporting the geometry defined, respectively as

    T_μν^(n)=σ l_μl_ν,
          T_μν^(m)=(ρ
    +p)(l_μn_ν+l_νn_μ)+pg_μν.

where σ, ρ and p are null radiation density, energy density
and pressure of the perfect fluid, respectively. In this regard, l_μ and n_μ
are linearly independent future pointing null vectors as

    l_μ=(1,0,0,0),     &   
    n_μ=(1/2(1-m(t,r)/r),-1,0,0
    ).

satisfying the following conditions

    l_μl^μ=n_μn^μ=0,    &    l_μn^μ=-1.

By this null vectors, the non-vanishing components of the total energy-momentum
tensor will be

    T_00   =   σ+ρ(1-m(t,r)/r),
    
    T_01   =   -ρ, 
    
     T_22   =   pr^2,
    
    T_33   =   pr^2 sin^2θ.

Moreover, using the metric ansatz (<ref>), we obtain

    𝒦^μ_ν=diag(0, 0,
    c/r,c/r).
 Consequently, we find

    (𝒦^2)^μ_ν=𝒦^μ_α𝒦^α_ν=diag(0, 0, c^2/r^2,c^2/r^2),
          (𝒦^3)^μ_ν=𝒦^μ_α𝒦^α_β𝒦^β_ν=diag(0, 0, c^3/r^3,c^3/r^3),
          (𝒦^4)^μ_ν=𝒦^μ_α𝒦^α_β𝒦^β_λ𝒦^λ_ν=diag(0, 0,
    c^4/r^4,c^4/r^4),

as well as the following required quantities

    [𝒦]   =   𝒦^μ_μ=2c/r,
    [𝒦^2]   =   (𝒦^2)^μ_μ=2c^2/r^2,
    [𝒦^3]   =   (𝒦^3)^μ_μ=2c^3/r^3,
    [𝒦^4]   =   (𝒦^4)^μ_μ=2c^4/r^4.

Now, using the equations (<ref>) and (<ref>), we can find the  𝒰_is in the equation (<ref>) as

    𝒰_1   =   2c/r, 
    𝒰_2   =   2c^2/r^2, 
    𝒰_3   =   0, 
    𝒰_4   =   0.

Using the equations (<ref>), (<ref>), (<ref>) and (<ref>), we can obtain the non-vanishing components of the massive gravity
term χ_μν in the field equation
(<ref>) as

    χ_00   =   [c_1c/r+ c_2c^2/r^2]( 1-m/r),
    χ_01   =   χ_10=-1/r(c_1c+c_2c^2/r), 
    χ_22   =   -c_1c r/2,
    χ_33   =   -c_1 c r sin^2θ/2.

 Then, for the 00 component of the field equation (<ref>), we have

    1/r^3[ r ṁ +r m^' -m  m^']
    =σ+ρ(1-m/r)-ℳ^2[c_1c/r
    + c_2c^2/r^2]( 1-m/r),

where dot and prime signs denote the derivative with respect to time and
radial coordinates, respectively. For the 01 and 10 component of the field equation, we have

    -m^'/r^2=-ρ+ℳ^2/r(c_1c+c_2c^2/r).

Finally, for the 22 and 33 component of the field equation, we obtain

    -1/2r
    m^''=pr^2+ℳ^2c_1c r/2.




§ DYNAMICS OF THE COLLAPSING SYSTEM

In this section, we discuss the  collapsing system dynamics, which is a time dependent system,
in massive gravity. This can be done by finding a solution for the field equations obtained  in the previous section.
We assume that the matter field follows the barotropic equation of state,  which is given by

    p=kρ,

where k is the barotropic parameter. Using the equations (<ref>),
(<ref>), (<ref>) and (<ref>) we get a solution for m as follows,

    m(t,r)=r^1-2k/1-2kf_1(t)+f_2(t)-1/2ℳ^2cc_1r^2-ℳ^2c^2c_2r,    
    k≠ 1/2

where f_1(t) and f_2(t) are arbitrary functions of t
given by f_1(t)=ρ(t,r)r^2(1+k), and


    σ(t,r)=r^-(1+2k)/1-2kḟ_̇1̇(̇ṫ)̇+ḟ_̇2̇(̇ṫ)̇/r^2

with k≠ 1/2,  where dot represents derivative with respect to t.

Therefore, the metric given in the equation (<ref>) can be written as
 the generalized Vaidya metric in massive gravity with the metric function
    f(t, r)=-1+r^-2k/1-2kf_1(t)+f_2(t)/r-1/2ℳ^2cc_1r-ℳ^2c^2c_2,    
    k≠ 1/2.

Now we will investigate the existence  of naked singularity (NS) in generalized Vaidya
spacetime. This will be done  using the outgoing radial null
geodesics, which will end up in the past at a singularity. So, the
geodesics will terminate  in the central physical singularity
located at r=0. It is possible for this singularity to be either
a naked singularity or a black hole (BH).
 Now for a   locally naked singularity, such null
geodesics exist. Furthermore, if the singularity is not a naked singularity, then this system forms a black hole.
Thus, by analyzing the   radial null geodesics that   emerge from the
singularity, we can understand the nature of such a singularity.

It may be noted that a  singularity can be formed by a catastrophic gravitational
collapse. In general, such a singularity can be either  a naked singularity  or a black hole. However, in  general relativity, such a singularity
formed from a gravitational collapse is always a black hole. This is because of the
 cosmic censorship in general relativity.  So,  in general relativity, the
 singularity is always covered by an event
horizon. However, this need not be the case for a more general theory. In fact,  it is possible  for
inhomogeneous dust cloud  to form  a naked singularity through
a collapse <cit.>. It may be noted that some interesting results have been obtained
for fluids whose equations of state is different from
a dust cloud <cit.>. So, it is possible to generalize the cosmic censorship in general relativity <cit.>.

Now, let us consider  R(t, r ) as the physical radius at time t of the
shell at r. Using the  scaling freedom in this system, we can write  R(0,r)=r at   the starting
time  t=0. So, different shells become singular at
different times for  the inhomogeneous
case. It is possible for  future directed radial null
geodesics to come out  of the singularity, with a well defined
tangent at the singularity. So,  dR/dr must tend to a finite limit, as the system  approaches
the past singularity along
these trajectories.

 As  the singularity is formed at  R(t_0,0)=0, the matter shells are crushed to zero radius at (t_0,r)=0.
This  singularity, which is formed at r=0,  is called as the  central singularity.
Now it is possible for  future directed
non-space like curves to have  their past end points
at this singularity. Such a singularity would then be called as a naked singularity.
So, for such a system,  the outgoing null geodesics
terminate in the past at the central singularity located at r=0. This occurs  at
t=t_0, and for this point, R(t_0, 0)=0. So,  along these geodesics, we obtain
R→ 0 as r→ 0 <cit.>.

We can write the  equation for the outgoing radial null geodesics using
 the equation (<ref>), and  setting ds^2=0 and
dΩ_2^2=0. Thus,  we obtain

    dt/dr=2/(1-m(t,r)/r).

It may be noted  that this system has a singularity at  r=0, t=0. Now, if the function X is given by  X=t/r, then
we can study the limiting behavior of X as we
approach the singularity located at r=0, t=0,  along the radial null
geodesic. This limiting value of X will be denoted  by X_0, and so we can write

    [ X_0; ;     ][ =lim   X;     t→ 0;     r→ 0 ][ =lim  t/r;      t→ 0;      r→ 0 ][ =lim  dt/dr;        t→ 0;        r→ 0 ][ =lim  2/(1-m(t,r)/r).;                  t→ 0;                  r→ 0 ]

Using the equations (<ref>) and (<ref>), we have

    2/X_0=
    [  lim; t→ 0; r→ 0 ][1-r^-2k/1-2kf_1(t)-f_2(t)/r+1/2ℳ^2cc_1r+ℳ^2c^2c_2],      k≠ 1/2.

Now, choosing f_1(t)=α t^2k and f_2(t)=β t, we
obtain an algebraic equation in terms of X_0 from equation
(<ref>) as

    α/1-2kX_0^1+2k+β
    X_0^2-(1+ℳ^2c^2c_2)X_0+2=0,     
    k≠ 1/2.

where α and β are constants.
It may be noted that naked singularity can also from in  general relativity,
as this equation would have positive roots
even for ℳ= 0.   Such  naked singularities in Vaidya spacetime
have been studied  in  general relativity <cit.>-<cit.>. However,
in this paper, we will analyze the effects of the mass term on the formation
of naked singularity, and observe what effect can such a mass term have on the
formation of naked singularities. As can be observed the above equation,  a non-zero mass term
ℳ= 0
will change the positive roots of this equation.
So, the addition of such a term will change the effect the behavior of this system.
To analyze such a behavior we first observe that
  the expression of
f_1(t),   can either be a constant or a non-linear
function of t, depending on the EoS, i.e. the cosmological era.
In the early universe (k≥ 0), we see that f_1(t) grows
with t, whereas in the late universe (k<0), f_1(t) decays
with time. This fact is demonstrated in Fig.<ref>(a). On
the other hand f_2(t) is a linear function of t. Nature of
f_2(t) is demonstrated in Fig.<ref>(b). It can be seen
that these choices of the arbitrary functions f_1(t) and
f_2(t) are somewhat self-similar in nature. The choice of
f_1(t) is driven by the presence of r^-2k in the second
term in equation (<ref>). Similarly the choice of f_2(t) is
based on the presence of the term 1/r in the third term in the
equation (<ref>). These self-similar choices follow from the
definition of X_0 given in equation (<ref>). We did not
consider non self-similar cases in order to avoid computational
difficulties. Choices other than self-similar ones will leave
residual r or t coordinates in the second and third terms of
equation (<ref>). In the limiting condition this will either
result in elimination of terms or creation of mathematically
undefined terms, both of which are undesirable. So this can be
considered a special class of solution given by the self-similar
choice of the functions f_1(t) and f_2(t).




A black hole will be formed if we obtain only non-positive
solution of this equation. However, if we obtain a positive real
root for this equation, then this system will be described by a
naked singularity. Here it is difficult to find exact solutions
for X_0 except for some particular values. This is because the
governing equation of the system is a very complicated one. These
exact solutions are given in the following Table 1. It is clear
from the table that certain conditions between the parameters are
required to be satisfied in order to make the solutions positive.







 k      Regime               Solution 1          Solution 2 

                                     



 0     Dust     U-α-√((U-α)^2-8β)/2β     U-α+√((U-α)^2-8β)/2β 

                                      



 1    stiff
 fluid     β/3α+2^1/3Y/3α(Z+√(4Y^3+Z^2))^1/3-(Z+√(4Y^3+Z^2))^1/3/3×
2^1/3α                   -  

                                      




-1/2    DE     U-√(U^2-2β(α+4))/2β     U+√(U^2-2β(α+4))/2β 








Table 1: Exact Values of X_0 for specific values of the
EoS parameter k obtained from Eq.(<ref>) 

where 

                     U=1+c^2c_2ℳ^2,


                     Y=3α+3c^2c_2ℳ^2α-β^2,


                     Z=-54α^2+9αβ+9c^2c_2ℳ^2αβ-2β^3.





We now analyze the results furnished in the Table 1.


Case 1: k=0


This corresponds to the pressureless dust regime of the universe.
For the solutions to be real and finite, we must have β≠
0, (U-α)^2≥ 8β.


For positivity of the first solution we have for
β>0,  U-α≥√((U-α)^2-8β)⟹β≥ 0. But since β≠ 0, we must restrict ourselves
to β>0, which is in agreement with our assumption. So the
first solution is positive for any positive β. For
β<0, we have for positive solution
U-α≤√((U-α)^2-8β)⟹β≤ 0. But since β≠ 0, we must have β<0
which is in agreement with our assumption. So, for β<0,
the solution 1 is always positive and represents a NS. Hence, for any
non-zero real β the first solution represents an NS.


For solution 2 to be positive we must have for β>0,
U-α≥
-√((U-α)^2-8β)⟹β≥
0. This gives β>0, just like the previous case. Hence, the
solution is positive. Similarly for β<0 case also we get
positive solution. Thus, this solution also represents NS.
Therefore for k=0, we get NS as the end state of collapse.




Case 2: k=1


This corresponds to the early stiff fluid era of our universe.
Here, the situation is much more chaotic mathematically. We get
only one solution, which turns out to be a relatively complicated
one. Now physically speaking, in the early universe, due to big
bang extreme amount of chaos is expected. Moreover, there are
quantum fluctuations, so mathematically the scenario is justified.
Here, in order to have a positive solution we should have
β/3α+2^1/3Y/3α(Z+√(4Y^3+Z^2))^1/3≥(Z+√(4Y^3+Z^2))^1/3/3×
2^1/3α. Moreover for negative Y, Z^2≥ 4Y^3 for
the solution to be real.




Case 3: k=-1/2


This represents the dark energy era corresponding to the late time
accelerated expanding universe. Here, the solution to be real and
finite we have, β≠ 0, U^2≥
2β(α+4). For positive solution, we must have
β(α+4)≥ 0 in which for β>0 ⟹
α≥ -4. For β<0, we should have
β(α+4)≤ 0 ⟹ α≤
-4 in order to get a positive solution.




               Numerical Solutions of X_0 and their interpretations


In order to understand the dynamics of collapse, we need to have a
knowledge of X_0 not at discrete points of k, but throughout
the cosmologically meaningful region k<1, i.e., from early to
late universe. To achieve this we proceed to obtain numerical
solutions of X_0, by assigning different initial conditions to
the parameters describing this system. To visualize these
solutions we obtain contours for k-X_0 for different numerical
values of the involved parameters.

It may be noted from  the plots that the trajectories run across
the positive range of X_0 thus confirming the formation of NS.
In Fig.<ref>(a), we can observe the dependence of X_0 on
the EoS parameter k for different values of the massive gravity
parameter ℳ. We see that an increase in the value of
ℳ decreases the tendency of formation of NS. Hence we
observe that the dynamics of the system gets deformed by the
addition of graviton mass to this system. In Fig.<ref>(b),
the k-X_0 trajectories for different values of α are
obtained. Here also different values of α deform the
dynamics of this system. Greater the value of α, greater is
the tendency to form NS.










In Fig.<ref>(a), we observe the effect of c on the
collapsing system. It is observed that greater the value of c,
lesser is the tendency to form NS. In Fig.<ref>(b), we can
observe the effect of c_2 on the system. Here also we see that
an increase in c_2 decreases the possibility of NS.

















In Figs.<ref>, <ref> and <ref>, 3D-plots are
obtained to get a more comprehensive view of the dynamics of the
collapse. In all the figures the resulting surfaces entirely lie
in the positive half-space of X_0, thus showing the presence
of NS. In Fig.<ref>(a), the variation of k-ℳ
surfaces are obtained against X_0. We see that in the dark
energy regime (k<-1/3) the surface pushes towards the positive
direction of X_0, accompanied with a decrease in
ℳ, thus showing an increased tendency to form NS. In
Fig.<ref>(b), k-α surfaces are obtained against
X_0. Here also in the dark energy regime, there is an
increased tendency of NS accompanied by an increase in α.
In Fig.<ref>(a), k-c surface is obtained against X_0.
Here accompanied by a decrease in the value of c, we witness an
increased tendency of NS in the dark energy regime. In
Fig.<ref>(b), k-c_2 surface is obtained against
X_0. The results obtained are same as that of
Fig.<ref>(a). In Fig.<ref>, k-β surface is
obtained against X_0. We see that the k-β surface is
parallel to the β axis. This shows that the system is not
deformed by β and hence the collapse dynamics does not
depend on it. This is an important result. Finally just like the
previous cases here also the surface gets lifted in the dark
energy regime towards the positive direction of X_0 axis.




 Strength of Singularity


 It is important to know about the destructive capacity of a singularity, and this is measured using the 
 concept of  strength of singularity. The strength of singularity is related to the extension 
 of spacetime through the singularity. Now this can be quantified  using the Tipler's  formalism <cit.>-<cit.>. 
 Now using the  Tipler's formalism  <cit.>-<cit.>, the 
condition for a singularity to be strong is given by,

    [ S=lim  τ^2ψ;        τ→ 0;             ][ =lim  τ^2R_μνK^μK^ν>0;                  τ→ 0;                       ]

where R_μν is the Ricci tensor. Here ψ is a scalar
given by ψ=R_μνK^μK^ν,  
K^μ=dx^μ/dτ is the tangent to the non-spacelike geodesics at the singularity,  
and τ is the affine
parameter.  It has been demonstrated  that <cit.>, 

    [ S=lim  τ^2ψ;        τ→ 0;             ][ =1/4X_0^2(2ṁ_̇0̇);                  ;                   ]

where m_0 is given by 

    [ m_0=lim   m(t,r);             t→ 0;            r→ 0. ]

Furthermore, it is also possible to write 

    [ ṁ_̇0̇=lim  ∂/∂ t(m(t,r));                     t→ 0;                    r→ 0. ]


Using Eq. (<ref>) in the above relation (<ref>), we obtain 

    [ S=lim  τ^2ψ;        τ→ 0;             ][ =1/4X_0^2[2kα/1-2kX_0^2k-1+β];                              ;                               ]

It may be noted that it has been demonstrated that  X_0 is related to
the limiting values of mass as <cit.>

    X_0=2/1-2m_0'-2ṁ_̇0̇X_0

where m_0' is given by 

    [ m_0'=lim  ∂/∂ r(m(t,r));                    t→ 0;                   r→ 0. ]

Here ṁ_̇0̇ is given by the Eq. (<ref>). Now using
Eqs. (<ref>), (<ref>) and (<ref>) in
Eq. (<ref>), we obtain 

    2α/1-2kX_0^2k+1+2β
    X_0^2-(1+2ℳ^2c^2c_2)X_0-2=0

In a particular case, if k=-1/2 (dark energy) is considered, we
obtain a solution for Eq. (<ref>) as

    X_0=0.25/β(1+2c^2c_2ℳ^2+√(8β(2-α)+(1+2c^2c_2ℳ^2)^2))

which is positive when  1+2c^2c_2ℳ^2≥
2√(2β(α-2)) and α>2 for β>0, 
provided   c_2>0. The existence of such
a positive root signifies that the singularity is naked. Using the
above value of X_0 in the Eq. (<ref>), we get

    [ S=lim  τ^2ψ;        τ→ 0;             ][ =0.03125/βu^2(1-8αβ/u^2);                         ;                          ]

where we have 
u=1+2c^2c_2ℳ^2+√(8β(2-α)+(1+2c^2c_2ℳ^2)^2).

Here, we can write 

    [ S=lim  τ^2ψ >0;           τ→ 0;                ]

for 8αβ < u^2 for β>0. This is the condition
for a strong naked singularity. A 3D plot for S is
shown in Fig.<ref> for a particular scenario. The plot shows
that the surface lies in the positive region thus giving a strong
naked singularity.













§ THERMODYNAMICS

In this section, we would like to study the thermodynamics of  generalized Vaidya spacetime in massive gravity.
The thermalization temperature, for such a spacetime, is given by the following relation <cit.>,

    T=1/4πd/drf(t,r)|_r=r_h,

where r_h is the event horizon obtained from the following relation (see the equation (<ref>)),

    -1+r^-2k/1-2kf_1(t)+f_2(t)/r-1/2ℳ^2cc_1r-ℳ^2c^2c_2=0.

Real positive root of the above equation gives the event horizon
radius. In Fig.<ref>, we can see the typical behavior of
f(t,r) in terms of r. We show that it is possible to have two radii at
which f(t,r)=0, and the bigger one (solid red) shows the event
horizon radius (about r=4 of Fig.<ref>). We can also see
that the increasing value of k increases the value of outer
event horizon radius. In the case of k>0.5 we can see only one zero (see red and blue lines). Also, we can see extremal case with k=-2
and ℳ=0.287. We should note that having one or two
horizons is a function of the massive parameter ℳ.






In the special case of k=-1, one can obtain real positive root as,

    r_h=1/2f_1(Y^1/3+c_1^2c^2ℳ^4+4f_1(c_2c^2ℳ^2+1)/Y^1/3-c_1cℳ^2),

where we defined,

    Y=2√(3)f_1√(M)-12f_2f_1^2-6( c_2c^2ℳ^2+1 ) ℳ^2c_1cf_1-ℳ^6c_1
    ^3c^3,

where

    M   =   -3c_1^2c_2^2c^6ℳ^8+ ( -16 f_1c_2^3c^6+6f_2c_1^3c^3-6c_1^2c_2c^4) ℳ^6
       +   (36 f_1f_2c_1c_2c^3-48 f_1c_2^2c^4-3c_1^2c^2) ℳ^4
       +   36 f_1c ( f_2c_1-4/3c_2c ) ℳ^2+36f_1^2f_2^2-16 f_1.

By using the equation (<ref>) one can obtain,

    T=1/4π(ℳ^2/2c_1c-2kf_1/1-2kr_h^-2k-1-f_2r_h^-2).

In the plots of the Fig.<ref> we can see typical behavior of the temperature for some values of k in terms of r (a) and in terms of ℳ (b). Fig.<ref> (a) plotted in terms of r, however from the Figs.<ref> we show that selected value of parameters yields r_h≈4. Instead of small radius, in the special case of k=-1, it is approximately linear function of radius. We can see that for the case of k<0.5, temperature is increasing function of r to yields a constant for the large radius. Situation is vise for the cases of k>0.5 (see red and blue solid lines).

Then, by using the relation (<ref>) one can obtain temperature in terms of ℳ in which we can find it as increasing function of ℳ. Also, we can see that increasing β (f_2) decreases the temperature.






Now, we can write entropy as,

    S=π^2 r_h^2,

where we used π G=1. Hence, we can use the following relation to calculate total energy,

    U=∫TdS,

which yields to the following expression,

    U=-1/2 π f_2ln( r_h) +1/8ℳ^2π c_1cr_h^2-π k f_1/4k^2-4k+1r_h^1-2k.

In the Fig.<ref> (a) we can see typical behavior of U for some values of k and find that value of k reduces value of the internal energy. Also, in the Fig.<ref> (b) we can see that internal energy is increasing function of mass parameter. Internal energy may be used to obtain Helmholtz free energy.





Helmholtz free energy can be obtained via the following relation,

    F=U-TS,

which yields to the following expression,

    F=kf_1r_h^1-2k+√(k-1/2)[kf_1r_h^1-2k+2f_2(k-1/2)(lnr_h-1/2)]/-4/π√(k-1/2).

It is interesting to note that Helmholtz free energy is independent of ℳ.
In the plots of the Fig.<ref>, we can see typical behavior of the Helmholtz free energy in terms of radius for various values of k.





Finally, we can study specific heat in constant volume,

    C=(dU/dT)_V,

which yields to the following expression,

    C=π r_h^2kf_1r_h^1-2k-(f_2-1/2ℳ^2c_1cr_h^2)(k-1/2)/f_2(k-1/2)-kf_1(k+1/2)r_h^1-2k,

In the plots of the Fig.<ref> we can see typical behavior of the specific heat in terms of the mass parameter and radius. We can see that specific heat may be positive or negative (for k of order unit) which means some instability with possible phase transition.  One can study such instability in the context of thermal fluctuations <cit.> and find that presence of thermal fluctuations may remove mentioned instabilities. Also, Fig.<ref> (b) shows that specific heat is increasing function of mass parameter.











§ CONCLUSIONS AND DISCUSSION

In this paper, we have analyzed the gravitational collapse in
massive gravity. It was observed that the dynamics of this system
are changed due to the addition of a mass to this system. In this
paper, we first obtained  equation of motion for a time dependent
solution  in  massive gravity. Then barotropic equation of state
was used to further analyze such solutions. Finally, we applied
this solution to a dynamics of gravitational collapse. It was
observed that the gravitational collapse depends on the value of
the mass used to deform this theory.

Contours for X_0 are obtained against the barometric parameter
k for different values of other parameters like ℳ,
α, β, etc. Various regimes of the fluid content of
the universe has been plotted such as radiation (k>0),
pressure-less dust (k=0), dark energy (k<0), phantom (k<-1).
From the figures we see that there is no trajectories in the
negative region. This rules out the existence of black holes as
the end state of collapse in the context of massive gravity with the considered
parameters.

From the first Figs.<ref>, <ref>, <ref>,
<ref>, <ref>, we see that the contours and surfaces
of X_0 obtained against the parameters  ℳ,
α, β and k lie in the positive region. This
confirms the existence of positive roots of equation (<ref>)
and indicates the end state of the collapse can be a  naked
singularity in the context of massive gravity with the considered
parameters.

In Fig.<ref>(a), we  see that an increase in the value of
massive gravity parameter ℳ decreases the tendency of
formation of NS. Hence, we observe that the dynamics of the system
gets deformed by the addition of graviton mass to this system. In
Fig.<ref>(b), the k-X_0 trajectories for different
values of α are obtained and represents that greater values
of α increase the tendency to form NS. In contrast, in
Figs.<ref>(a) and <ref>(b), we observe that greater
values of c and c_2 decrease the tendency to form the NSs.
In Fig.<ref>(a), the variation of k-ℳ surfaces
are obtained against X_0. We see that in the dark energy
regime (k<-1/3) the surface pushes towards the positive
direction of X_0, accompanied with a decrease in
ℳ, thus showing an increased tendency to form NS. In
Fig.<ref>(b), k-α surfaces are obtained against
X_0. Here also in the dark energy regime, there is an
increased tendency of NS accompanied by an increase in α.
In Fig.<ref>(a), k-c surface is obtained against X_0.
Here accompanied by a decrease in the value of c, we witness an
increased tendency of NS in the dark energy regime. In
Fig.<ref>(b), k-c_2 surface is obtained against
X_0. The results obtained are same as that of
Fig.<ref>(a). In Fig.<ref>, we see that the k-β
surface is parallel to the β axis around k=0 and this
shows that the system is not deformed by β and hence the
collapse dynamics does not depend on it in the early universe
regime. But what matters is that the surface remains totally in
the X_0>0 region. Eventually the surface pushes up towards the
positive axis for k≤ -0.5, i.e. in the dark energy regime.
This shows an increased tendency to form naked singularities. In
the figs. 4-6, we see that the surfaces totally lie in the
X_0>0 region thus showing the possibility of formation of
naked singularities. We have also studied the strength of the
singularity formed and showed that under certain conditions, and concluded that 
strong singularities are possible in massive gravity. This is illustrated in Fig.<ref>.
It would be interesting to analyze the strength of such  naked singularities, 
for various models, and compare the results of massive gravity with general relativity. 

From the above discussion, we conclude that this work provides
counterexamples of the cosmic censorship conjuncture, which states
that every singularity must be covered by an event horizon, in the
context of massive gravity.  We would like to point out that there
are two forms of the cosmic censorship conjuncture
<cit.>-<cit.>. According to the  strong cosmic
censorship conjuncture, no locally naked singularities can occur.
However, according to the weak cosmic censorship conjuncture,
singularities can be locally naked, but they cannot be  globally
naked. It is also possible to analyze the strength of singularities, and this can be done 
using the Tipler's formalism   <cit.>-<cit.>. 
We have applied this formalism to the massive gravity, and demonstrated 
that it possible to have a strong naked singularity 
in massive gravity. 


Finally, we study thermodynamics of the model and calculate some
thermodynamical quantities to investigate effect of mass
parameter. For example, in Figs.<ref>(b) and <ref>(b),
we find that the thermalization temperature and specific heat
respectively are  increasing function of ℳ. We also
found some instabilities, corresponding to special values of k,
and we found stable/unstable black hole phase transition. For the
future work we would like to focus on the instabilities and
consider effect of thermal fluctuations to see that what can
happen with the instable regions.





§ ACKNOWLEDGMENTS


P. Rudra acknowledges University Grants Commission (UGC),
Government of India for providing research project grant (No.
F.PSW-061/15-16 (ERO)). P. Rudra also acknowledges Inter
University Centre for Astronomy and Astrophysics (IUCAA), Pune,
India, for awarding Visiting Associateship. F. Darabi acknowledges
financial support of Azarbaijan Shahid Madani University (No.
S/5749-ASMU) for the Sabbatical Leave, and thanks the hospitality
of ICTP (Trieste) for providing support during the Sabbatical
Leave. Y. Heydarzade acknowledges the support of Azarbaijan Shahid Madani University under the approvement (No. 214/D/24019-ASMU). 







99


superA.G. Riess et al., Astron. J. 116, 1009 (1998)
super1S. Perlmutter et al., Nature 391 , 51 (1998)
super2A. G. Riess et al., Astron. J. 118, 2668 (1999)
super0S. Perlmutter et al., Astrophys. J. 517, 565 (1999)
super4A. G. Riess et al., Astrophys. J. 560, 49 (2001)
super5 J. L. Tonry et al., Astrophys. J. 594, 1 (2003)
DEP. J. E. Peebles, and B. Ratra, Rev. Mod. Phys.75, 559 (2003)
DE1E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys.D15, 1753 (2006)
DE2J. A. Frieman, M. S. Turner, and D. Huterer,  Annual Review of Astronomy and Astrophysics, 46, 385 (2008)
M1
M. Khurshudyan, B. Pourhassan, A. Pasqua, Can. J. Phys. 93 (2015) 449
2a
H. van Dam and M. J. G. Veltman, Nucl. Phys. B 22, 397 (1970)

a2
Y. Iwasaki, Phys. Rev. D 2, 2255 (1970)

M2
S. Upadhyay, B. Pourhassan, H. Farahani, Phys. Rev. D 95, 106014 (2017)
1  W. Pauli  and M.  Fierz Helv. Phys. Acta 12, 297 (1939)
2
M. Fierz Helv. Phys. Acta 12, 3 (1939)



4 A. I. Vainshtein, Phys. Lett. B
39, 393
(1972)
5 E. Babichev and C. De ayet, Class. Quant. Grav.
30, 184001
(2013)
6 D. G. Boulware, S. Deser, Phys. Rev. D
6, 3368
(1972)
7de  Rham  C,  Gabadadze  G  and  Tolley  A  J  2011
Phys. Rev. Lett.
106,
231101 (2011)
8 de Rham C and Gabadadze G
Phys. Rev.
D
82,
04402 (2010)
9  de Rham C, Gabadadze G and Tolley A J
Phys. Lett.
B
711,
190 (2012)
10  S. F. Hassan, R. A. Rosen and A. Schmidt-May, JHEP
1202, 026
(2012)
11S. F. Hassan, A. Schmidt-May and M. von Strauss, Phys. Lett. B
715, 335
(2012)
12S. F.  Hassan  and R. A. Rosen
Phys. Rev. Lett.
108,
041101 (2012)
13 S. F.  Hassan S F and R. A. Rosen  JHEP
1204,
123 (2012)
14  K. Hinterbichler, Rev. Mod. Phys.
84, 671
(2012)
lore
  A. H. Chamseddine and V. Mukhanov,  JHEP  1208, 036 (2012)
  l  I. Arraut,
   arXiv:1505.06215 [gr-qc].
  l1S. Dengiz,  arXiv:1409.5371 [hep-th]
lore2 G. Goon, K. Hinterbichler, A. Joyce and M. Trodden,
   JHEP   1507, 101 (2015)


   gaus S. H. Hendi, S. Panahiyan, B. E. Panah,   JHEP 01, 129 (2016)
   gauss Se. H. Hendi, G. Q. Li, J. X. Mo, S. Panahiyan and B. E. Panah,
    Eur. Phys. J. C 76, 571 (2016)
    phase S. H. Hendi, S. Panahiyan, B. E. Panah and M. Momennia,
    Ann. der Phys.  528, 819 (2016)

cosm12  M. Wyman, W. Hu and P. Gratia,   Phys. Rev. D 87, 084046 (2013)
cosm14 M. S. Volkov,  Phys.Rev. D 90, 024028 (2014)

ads
A. Sinha,   JHEP 1006, 061 (2010)
M3
J. Sadeghi and B. Pourhassan, JHEP12 (2008) 026
M4
B. Pourhassan  and J. Sadeghi, Can J Phys 91 (2013) 995
M5
J. Sadeghi, B. Pourhassan and S. Heshmatian, Advances in High Energy Physics 2013 (2013) 759804

cft
V. Niarchos,  Fortsch. Phys. 57, 646 (2009)

ee X. X. Zeng, H. Zhang and  L. F. Li, Phys. Lett. B 756,  170 (2016)
cc W. J.  Pan and Y. C. Huang,   arXiv:1612.03627
conduct  L. Alberte and A. Khmelnitsky, Phys. Rev. D 91, 046006 (2015)
time V. Ziogas, JHEP 09, 114 (2015)
time1  V. Keranen and P.  Kleinert, Phys. Rev. D 94, 026010 (2016)
sphe P. C. Vaidya, Curr. Sci. 12, 183 (1943)
sphe1 P. C. Vaidya, Nature 171, 260 (1953)
td12 P. Li, X. Z. Li and X. H. Zhai, Phys. Rev. D 94, 124022 (2016)
td21
Y. P.  Hu, X. X.  Zeng and  H. Q. Zhang, Phys. Lett. B 765,  120 (2017)
M6
K. B. Fadafan, B. Pourhassan and J. Sadeghi, Eur. Phys. J. C 71, 1785 (2011)
Elena
E. Caceres, A. Kundu and D. L. Yang, JHEP 1403, 073 (2014)
gc M. Sharif and A.  Siddiqa,  Gen. Rel. Grav. 43, 73 (2011)
gc12  M. Sharif and A. Siddiqa,  Mod. Phys. Lett. A 25, 2831 (2010)
r1 P. Rudra, R. Biswas, U. Debnath,  Astrophys. Space Sci. 335, 505 (2011)
r2 U. Debnath, P. Rudra, R. Biswas,  Astrophys. Space Sci. 339, 135 (2012)
r3 P. Rudra, R. Biswas, U. Debnath,  Astrophys. Space Sci. 354, 597 (2014)
r4 P. Rudra, U. Debnath, Can. J. Phys. 92(11), 1474 (2014)
r5 P. Rudra, M. Faizal, A. F. Ali,  Nucl. Phys. B. 909, 725 (2016)
r6 Y. Heydarzade, P. Rudra, F. Darabi, A. F. Ali, M.
Faizal,  Phys. Lett. B. 774, 46 (2017)

refe1 I.  H.  Dwivedi and P.  S.  Joshi, Class. Quant. Grav. 6, 1599 (1989)
referen  I.  H.  Dwivedi and P.  S.  Joshi, Class. Quant. Grav. 8, 1339 (1991)
referen1 I.  H.  Dwivedi and P.  S.  Joshi, J. Math. Phys. 32, 2167 (1991)
referen2 I.  H.  Dwivedi and P.  S.  Joshi, Phys. Rev. D 45, 2147 (1992)
reference4 K. Lake, Phys. Rev. D 43, 1416 (1991)
reference4 S.  M.  Wagh and S.  D.  Maharaj, Gen. Rel. Grav. 31, 975 (1999)
reference4   S.  G.  Ghosh and A.  Beesham, Phys. Rev. D 61, 067502 (2000)
reference4   S.  G.  Ghosh and N.  Dadhich,, Phys. Rev. D 64, 047501 (2001)
refe2A.  Beesham and S.  G.  Ghosh,  Int. J. Mod. Phys. D12, 801 (2003)




reference R. G. Cai, Y. P. Hu, Q. Y. Pan and Y. L. Zhang, Phys. Rev. D 91 024032 (2015)
Eardley1 D. M. Eardley,  and  L. Smar ,  Phys. Rev. D 19,  2239 (1979)
Joshi1 P. S. Joshi, I. H. Dwivedi,  Commun. Math. Phys.  146,  333 (1992)
Joshi2 P. S. Joshi, T. P. Singh,  Phys. Rev. D 51,  6778 (1995)
Singh1 T. P. Singh, P. S. Joshi,  Class. Quant. Grav. 13, 559 (1996)

Tipler F. J. Tipler, Phys. Lett. A. 64, 8 (1977)
conj15V. D. Vertogradov, J  Phys. Conf. Series,  769, 012013 (2016)
conj16 V.D. Vertogradov, Grav. Cosmol. 22, 220 (2016)
 Maharaj M. D. Mkenyeleye, R. Goswami, and S. D. Maharaj, Phys. Rev. D 90,  064034 (2014)



T1
B. Pourhassan and M. Faizal,  Nucl. Phys. B 913, 834 (2016)
T2
J. Sadeghi, B. Pourhassan and M. Rostami, Phys. Rev. D 94, 064006 (2016)
T3
B. Pourhassan and M. Faizal, Phys. Lett. B 755, 444 (2016)
T4
B. Pourhassan, M. Faizal and U. Debnath, Eur. Phys. J. C 76, 145 (2016)
T5
M. Faizal and B. Pourhassan, Phys. Lett. B 751, 487 (2015)
T6
B. Pourhassan and M. Faizal, Europhys. Lett. 111, 40006 (2015)

conj12 T. P. Singh,   gr-qc/9606016
conj14 S. S. Deshingkar, S. Jhingan and  P. S. Joshi, Gen. Rel. Grav. 30, 1477 (1998)





















