Thermodynamics and reentrant phase transition for logarithmic nonlinear charged black holes in massive gravity
S. Rajaee Chaloshtary M. Kord Zangeneh S. Hajkhalili A. Sheykhi S. M. Zebarjad




Physics Department and Biruni Observatory, Shiraz University, Shiraz 71454,
Iran

mkzangeneh@scu.ac.ir
Physics Department, Faculty of Science, Shahid Chamran University of Ahvaz,
Ahvaz 61357-43135, Iran

Physics Department and Biruni Observatory, Shiraz University, Shiraz 71454,
Iran

asheykhi@shirazu.ac.ir
Physics Department and Biruni Observatory, Shiraz University, Shiraz 71454,
Iran
Max-Planck-Institute for Gravitational Physics (Albert-Einstein-Institute),
14476 Potsdam, Germany

zebarjad@shirazu.ac.ir
Physics Department and Biruni Observatory, Shiraz University, Shiraz 71454,
Iran


We investigate a new class of (n+1)-dimensional topological black hole
solutions in the context of massive gravity and in the presence of
logarithmic nonlinear electrodynamics. Exploring higher dimensional
solutions in massive gravity coupled to nonlinear electrodynamics is
motivated by holographic hypothesis as well as string theory. We first
construct exact solutions of the field equations and then explore the
behavior of the metric functions for different values of the model
parameters. We observe that our black holes admit the multi-horizons caused
by a quantum effect called anti-evaporation. Next, by calculating the
conserved and thermodynamic quantities, we obtain a generalized Smarr
formula. We find that the first law of black holes thermodynamics is
satisfied on the black hole horizon. We study thermal stability of the
obtained solutions in both canonical and grand canonical ensembles. We
reveal that depending on the model parameters, our solutions exhibit a rich
variety of phase structures. Finally, we explore, for the first time without
extending thermodynamics phase space, the critical behavior and reentrant
phase transition for black hole solutions in massive gravity theory. We
realize that there is a zeroth order phase transition for a specified range
of charge value and the system experiences a large/small/large reentrant
phase transition due to the presence of nonlinear electrodynamics.


04.70.-s, 04.30.-w
Thermodynamics and reentrant phase transition

for logarithmic nonlinear charged black holes in massive gravity
    S. M. Zebarjad
    December 30, 2023
===============================================================================================================





§ INTRODUCTION


The initial version of General Relativity (GR) presented by Einstein is a
theory which predicts the existence of a massless spin-2 particle called
graviton <cit.>. One of the most natural generalizations of massless GR
theory is to admit the existence of the massive graviton and hence
constructing a massive theory of gravity. According to the particle physics,
Higgs mechanism bestows mass to carrier particles of the electroweak
interaction. Strong motivation for studying the massive gravity models comes
from recent developments in gravitational waves observatory. The recent
observation of the the gravitational wave event GW170817 and of its
electromagnetic counterpart GRB170817A, from a binary neutron star merger,
confirmed that the speed of gravitational waves deviates from the speed of
light by less than one part in 10^15 <cit.>. This implies that
one should consider an upper bound for the mass of graviton <cit.>.
Thus the observation of the gravitational waves put severe constraints on
several theories of modified gravity. Any modified gravity model predicting
the speed of gravitational waves less than the speed of light must now be
seriously reconsidered <cit.>. From holographic point of view, massive
gravity could make momentum dissipative in dual systems <cit.>. Giving
mass to graviton state breaks the diffeomorphism invariance of stress-energy
tensor of dual theory holographically which in turn cause momentum
dissipation <cit.>.

Historically, the first massive gravity theory was derived by Fierz and
Pauli by adding the interaction terms to the linearized level of GR <cit.>. Their theory suffered a discontinuity in the predictions. For
instance, Van Dam, Veltman and Zakharov found a discontinuity in the
Newtonian potential in massless limit, recognized as VDVZ problem <cit.>
. Boulware and Deser suggested a nonlinear model <cit.>, while
Vainshtein introduced the linearized theory as culprit for VDVZ problem <cit.>. However, their theory was plagued with a ghost appeared as a 6th
degree of freedom. A four-dimensional covariant nonlinear theories of
massive gravity which are ghost-free in the decoupling limit to all orders
was established by de Rham-Gabadadze-Tolley (dRGT) <cit.>. Their
theories resum explicitly all the nonlinear terms of an effective field
theory of massive gravity <cit.>. Actually, the absence of higher
derivative terms in equations of motion does not allow the ghost to exist.
This theory has been explored from different point of view and known as the
only candidate for pathology-free nonlinear massive gravity so far <cit.>. Although, finding exact solutions is not easy in this model due to
the complexity of equations, it has always been an interesting and well
motivated topic of research in theoretical physics <cit.>.

Black hole solutions in the context of massive gravity have been explored
from different points of view. A Schwarzschild-de Sitter-like vacuum
solution of dRGT model has been found in <cit.>. A nontrivial black
hole solution with a negative cosmological constant has been obtained in 
<cit.>. It has been shown that the mass of the graviton in massive
theory is equivalent to the lattice in the holographic conductor model <cit.> (see also <cit.>). Charged black hole solutions of massive
gravity theory with negative cosmological constant have been constructed in 
<cit.> and their thermodynamic properties and phase structure have
been studied in both canonical and grand canonical ensembles. Exact
three-dimensional asymptotically AdS-like solutions in massive gravity have
been explored in <cit.>. The authors of Ref. <cit.> found the
classical solutions for stars and other compact objects in massive gravity
following the Vainshtein mechanism. Neutron stars in the context of this
theory have been investigated in <cit.>.

The importance of massive gravity does not restrict to the non-perturbative
aspects of gravity or black hole solutions. In the cosmological framework,
the massive gravity can be regarded as an alternative for the late time
cosmic acceleration <cit.>. Thanks to Vainshtein mechanism,
massive gravity is equivalent to GR at small scales, although in large
distances (i.e. infrared regime) it amends gravity. Thus, the nature of the
cosmological constant may change since it is the most universal infrared
possible source <cit.>. Furthermore, massive particles generate Yukawa
potential ∼e^-m r/r which has an exponential suppression
kicking in at the length scale ∼1/m. If we imagine that m∼ H
where H is Hubble scale, then the force of massive graviton would be
weakened at large scales and it may offer a nice solution for acceleration
problem <cit.>. More researches in this regard can be found in <cit.>.

When the linear Maxwell theory failed to explain the central singularity of
a point charge, the striking alternative was nonlinear electrodynamics
(NED). In 1930, Born and Infeld (BI) made the first efforts to construct a
NED theory <cit.>. Some years later, it turned out that their theory
plays a prominent role in D-brane physics. In the framework of particle
physics, NED was propounded by Heisenberg and Euler <cit.>. Its
extension in GR was done by Plebanski <cit.>. The theories of NED have
arisen a lot of attention after development in the string theory <cit.>
. For instance, loop calculations in open superstring theory lead to a low
energy effective action included BI type term <cit.>. It is noteworthy
to mention that at the high energy levels, due to the interactions with
other physical fields, real electromagnetic field cannot obey linear law.
Basically, it was argued that NEDs are the simplified phenomenological
descriptions of the pointed out interactions. In addition, NED can explain
the Rindler acceleration as a nonlinear electromagnetic effect <cit.>.
From quantum gravity viewpoint, NEDs are corrections to Maxwell field. So, a
large class of regular black holes (i.e. black holes without singularity)
have been studied in the literature. The pioneer effort in this direction is
attributed to Bardeen <cit.>.

The Lagrangian of a general theory of NED is assumed to be a function of the
Maxwell invariant ℱ=F_μνF^μν where F_μν
is electromagnetic tensor and not to contain any higher derivative terms of ℱ. In the weak field limit, these theories lead to linear Maxwell
one and so their expansions in terms of ℱ have the form ℒ(ℱ)=ℱ+O(ℱ^2). There has been
some efforts to construct BI type NED Lagrangian in literature. Two famous
ones are so-called Logarithmic (LNE) and Exponential (ENE) nonlinear
electrodynamics <cit.>. LNE removes divergences in the electric
field (like BI model), whereas the ENE just makes it weaker than Maxwell
theory. In the present work, we are interested in considering LNE
Lagrangian. It is notable to mention that logarithmic term of the field
strength may come out as an exact one-loop correction to the vacuum
polarization. This fact has been proved by Euler and Heisenberg when they
studied electrons in a background setup by a uniform electromagnetic field 
<cit.>.

LNE Lagrangian has been studied in various contexts. In Einstein-dilaton
theory, black hole solutions have been constructed and their thermodynamic
stability have been discussed <cit.>. AdS-dilaton black holes in the
presence of LNE have been investigated in <cit.>. Lifshitz-dilaton
black holes/branes coupled to LNE have been investigated in <cit.>.
Magnetic rotating dilaton strings and charged rotating dilaton black strings
in Einstein gravity with LNE electrodynamics have been found in <cit.>
and <cit.>, respectively. In the context of Lovelock gravity, authors
of <cit.> explored magnetic brane solutions in the presence of LNE.
Other studies on the black hole solutions in the presence of LNE can be
carried out in <cit.>.

The layout of this paper is as follows. In the next section, we introduce
the Lagrangian of massive gravity in the presence of cosmological constant
and LNE and construct exact black hole solutions. We calculate conserved and
thermodynamic quantities of the solutions in section <ref>. Thermal
stability of the solutions is checked in <ref>. In <ref>
, the critical behavior of (3+1)-dimensional black holes is examined. We
devoted the last section to summary and some closing remarks.



§ ACTION AND MASSIVE GRAVITY SOLUTIONS


The (n+1)-dimensional action of Einstein massive gravity with a negative
cosmological constant and NLE is <cit.>


    𝒮=1/16π∫ d^n+1x√(-g)[ ℛ
    -2Λ -8β ^2ln( 1+ℱ/8β ^2)
    +m^2∑_i^4c_i𝒰_i(g,Γ )] ,

where ℛ is the Ricci scalar, Λ =-n(n-1)/2l^2 with l
is the radius of AdS spacetime and ℱ=F_μνF^μν is
the Maxwell invariant and F_μν=2∂ _μA_ν ]
is the electromagnetic field tensor, while A_ν is the
electromagnetic vector potential. The constant β is the parameter of
nonlinearity so that the logarithmic nonlinear electrodynamics recovers the
linear Maxwell theory when β→∞. m^2 is the
positive massive gravity parameter so that the translational invariance is
recovered as m approachs to zero (m→ 0).

In action (<ref>), c_i's are constants and 𝒰_i's
are symmetric polynomials eigenvalues of (n+1)× (n+1) matrix 𝒦_ν^μ=√(g^μλΓ _λν) so
that


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

Here, rectangular braket [ 𝒦] =𝒦_μ^μ is the trace of the matrix square root which is defined as (√(
𝒦_ν^μ))× (√(𝒦_λ^v))=
𝒦_λ^μ. The dynamical metric is denoted by g and Γ is reference metric which is a 2-rank symmetric tensor.

One can obtain the equations of motion by varying the action (<ref>)
with respect to the metric tensor g_μν and gauge field A_μ
as

    G_μν+Λ g_μν+2F_μλF_ν^ 
    λ/1+ℱ/8β ^2+4β ^2ln( 1+
    ℱ/8β ^2) g_μν+m^2χ _μν=0,


    ∇ _μ[ ( 1+ℱ/8β ^2)
    ^-1F^μν] =0,

where, G_μν is the Einstein tensor and 

    χ _μν   =   -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).

In order to attain the static charged black hole solution, we consider the
line element of (n+1)-dimensional spacetime as 

    ds^2=-f(r)dt^2+f^-1(r)dr^2+r^2h_ij(x)dx^idx^j,             (i,j=1,2,...n-1)

where the function f(r) should be determined and h_ij(x)dx^idx^j
exhibits a hypersurface with constant curvature (n-1)(n-2)k where k
explicate the topology of event horizon or the boundary of t=constant - r=constant. It will be k=0 (R^n-1), k=1 (S^n-1) or k=-1 (H^n-1) . By taking k=0,1,-1, the black hole horizon can be flat
(zero), spherical (positive) or hyperbolic (negative) constant curvature
hypersurfaces with volume V_n-1. The reference metric can be considered
as <cit.> 

    Γ _μν=diag(0,0,c_0^2h_ij(x)),

where c_0 is a positive constant. By employing metric (<ref>) and (
<ref>) we can calculate 𝒰_i 's


    𝒰_1   =   (n-1)c_0/r,  
    𝒰_2   =   (n-1)(n-2)c_0^2/r^2,  
    𝒰_3   =   (n-1)(n-2)(n-3)c_0^3/r^3,  
    𝒰_4   =   (n-1)(n-2)(n-3)(n-3)c_0^4/r^4,

It is remarkable to note that 𝒰_3 and 𝒰_4 become
zero for (3+1)-dimensional spacetime and 𝒰_4 vanishes for (4+1)-dimensional spacetime. By considering the metric (<ref>), one
can integrate (<ref>) to calculate the electromagnetic field


    F_tr=-F_rt=2qr^1-n/1+Υ,

with 

    Υ =√(1+( q/β r^n-1) ^2).

In latter equations, q is a constant which has relation with total charge
of black hole. Substituting the reference metric (<ref>), (<ref>)
and electromagnetic field (<ref>) into the field Eqs. (<ref>), we
find


    f^'+(n-2)(f-k)/r+2Λ r/n-1-8β ^2r/
    n-1[ ln( 1+Υ) +1-Υ]
    -c_0m^2( c_1+(n-2)c_0c_2/r.       
    . +(n-2)(n-3)c_0^2c_3/r^2+
    (n-2)(n-3)(n-4)c_0^3c_4/r^3)    =   0,



    f^''+2(n-2)/rf^'+(n-2)(n-3)(f+k)/r^2
    +2Λ -8β ^2ln( 1+Υ/2) -(
    n-2) c_0m^2/r( c_1+(n-3)c_0c_2/r.       
    . +(n-3)(n-4)c_0^2c_3/r^2+
    (n-3)(n-4)(n-5)c_0^3c_4/r^3)    =   0,

where the prime denotes derivative with respect to r. Finally, one could
solve the above differential equations to obtain the metric function f(r)
as 

    f(r)    =   k-m_0/r^n-2-2r^2Λ/n(n-1)+8β
    ^2r^2/n( n-1) [ ln( 1+Υ/2
    ) +( 2n-1) ( 1-Υ) /n] 
    
          +8( n-1) q^2r^4-2n/n^2( n-2) 𝐅
    ( 1/2,-2+n/2(-1+n),4-3n/2-2n,
    -q^2r^2-2n/β ^2)  
          +c_0m^2r/n-1( c_1+(n-1)c_0c_2/r+
    (n-1)(n-2)c_0^2c_3/r^2+(n-1)(n-2)(n-3)c_0^3c_4/r^3
    ),

where 𝐅 is the hypergeometric function. Linear limit (β→∞) of the metric function is 

    f(r)    =   k-m_0/r^n-2-2Λ r^2/n(n-1)+2q^2/
    (n-1)(n-2)r^2(n-2)-q^4r^6-4n/4( n-1) (
    3n-4) β ^2
          +c_0m^2r/n-1( c_1+(n-1)c_0c_2/r+
    (n-1)(n-2)c_0^2c_3/r^2+(n-1)(n-2)(n-3)c_0^3c_4/r^3
    ) +O( 1/β ^4)

which is in accordance with the metric function of charged black hole in
Maxwell theory <cit.> along with a correction term corresponding to
the nonlinear electrodynamics. 


We also choose m_0 as a constant of integration which is related to
total mass of black hole. It is easy to obtain the parameter m_0 in (<ref>) by using Eq. f(r_+)=0. We find 

    m_0   =   kr_+^n-2-2r_+^nΛ/n(n-1)+8β
    ^2r_+^n/n( n-1) [ ln( 1+Υ _+/2
    ) +( 2n-1) ( 1-Υ _+) /n]
    
          +8( n-1) q^2r_+^-2n/n^2( n-2) 
    𝐅( 1/2,2-n/2(-1+n),4-3n/2-2n,
    -q^2r_+^2-2n/β ^2)  
          +c_0m^2r_+^n-1/n-1( c_1+(n-1)c_0c_2/r_+
    +(n-1)(n-2)c_0^2c_3/r_+^2+
    (n-1)(n-2)(n-3)c_0^3c_4/r_+^3) ,

in which Υ _+=Υ (r_+) and r_+ is the radius of
horizon which is the largest root of f(r_+)=0. As it is clear from Eq. (
<ref>), the position of r_+ completely depends on the metric
parameters. This fact is shown in Fig. <ref>. In this figure, we set l=k=1. Depending on the number of metric function's root(s), our solution
may be a black hole with two inner and outer horizons, an extreme black hole
or a naked singularity. As one can see from this figure, the number of
horizons, decreases by increasing the massive gravity parameter m and
electric charge q. It is notable to mention that a suitable choice of
metric parameters may result a Schwarzschild-like black hole (Fig. <ref>
). From Fig. <ref>, we observe that in the presence of massive gravity
and nonlinear electrodynamics, the mass parameter m_0 has non zero value
when r_+→ 0. Fig. <ref> also confirms the result of Fig. 
<ref>. Comparing dash-dotted and dashed as well as solid and
dash-double-dotted curves in Fig. <ref> shows that by increasing q
when m is fix, we could move from Schwarzschild-like black holes with just
one horizon to black holes with two horizons. One could also depict curves
to show the same behavior for m (as shown before in Fig. <ref>),
however we avoid this for economic reasons.

Fig. <ref> represents a significant behavior of the metric function f(r); the existence of several roots! The green (dashed) curve shows two
roots where the smaller one is extreme and the blue (dash-dotted) curve
exhibits the existence of three roots. Also, as one decreases the c_0
parameter, the largest root moves to larger values. Similar result has been
reported in <cit.>. It has been discussed that the physical reason
behind this property is anti-evaporation <cit.>. Anti-evaporation is a
quantum effect which causes the size of black hole to increase.

The gauge potential A_t can be calculated as 

    A_t( r)     =∫ F_rtdr  
        =μ +2r^nβ ^2/nq( Υ -1) +2(
    1-n) q^2r^2-n/n( n-2) 𝐅(1/2,
    n-2/2(n-1),4-3n/2(1-n),1-Υ ^2).

From holographic point of view, μ as a constant of integration is
chemical potential of quantum field that locates on the boundary. It can be
found by demanding the regularity condition on the horizon i.e. A_t(r_+)=0


    μ =2r_+^nβ ^2/nq( 1-Υ _+) -
    2( 1-n) q^2r_+^2-n/n( n-2) 𝐅(
    1/2,n-2/2(n-1),4-3n/2(1-n),1-Υ
    _+^2).

We wind up this section by calculating the Hawking temperature on the event
horizon


    T    =   f^'(r_+)/4π
       =   (n-2)k/4π r_++m^2c_0/4π r_+^3(
    c_0^3c_4(n-4)(n-3)(n-2)+c_0^2c_3(n-3)(n-2)r_++c_0c_2(n-2)r_+^2+c_1r_+^3)
    
          -r_+Λ/2π (n-1)+2r_+β ^2[ (
    1-Υ _+) +ln( 1/2( 1+Υ _+)
    ) ] /π (n-1).

In the next section, we will study the thermodynamics of our black hole
solutions.



§ THERMODYNAMICS OF MASSIVE GRAVITY SOLUTIONS 


In order to investigate the thermodynamics of our black hole solutions, we
have to calculate some thermodynamical quantities. We first obtain the
entropy of the black holes. Based on <cit.>, this thermodynamic
quantity is equal to one-quarter of the horizon area 

    S=r_+^n-1/4.

The Gauss law allows us to compute electric charge per unit volume V_n-1
using electromagnetic flux at infinity,

    Q= 1/4π∫ r^n-1( 1+ℱ/8β ^2
    ) ^-1F_μνn^μu^νdr,

in which the unit spacelike and timelike normals to the hypersurface of
radius r are respectively n^μ=( √(-g_tt))
^-1dt=( √(f( r) )) ^-1dt and u^ν=( 
√(g_rr)) ^-1dr=√(f(r))dr. Finally, the electric charge
per unit volume V_n-1 of black hole could be obtained as

    Q= q/4π.

Electric potential can measure at infinity with respect to the horizon using
the following definition 

    U=A_μχ ^μ| _r→∞-A_μχ ^μ| _r=r_+,

where χ =∂ _t is the null generator of the horizon. So it is
easy to show that electric potential is equal to chemical potential 

    U=μ .

The mass of our charged solution using the Hamiltonian approach can be found
as

    M=(n-1)m_0/16π.




Using (<ref>), (<ref>), (<ref>) and (<ref>), it is a
matter of calculation to show that the mass function is


    M(S,Q)=    -Λ (4S)^n/n-1/8nπ+m^2c_0/
    16π[ (n-1)(n-2)(4S)^n-3/n-1( c_3+(4S)^-1/
    n-1(n-3)c_4) c_0^3+(4S)^n-2/n-1
    (n-1)c_2c_0+4Sc_1]  
        +β ^2/2nπln[ 1/2+1/2√(1+
    Q^2π ^2/β ^2S^2 )] (4S)^n/n-1+
    8(n-1)^2 Q^2π/n^2(n-2)𝐅( 1/2,n-2/
    2n-2,3n-4/2-2n,-Q^2π ^2/β ^2S^2) (4S)^
    2-n/n-1
        +1/16n^2π[ 8β ^2(2n-1)+k n^2(n-1)(4S)^-2
    /n-1-8β ^2(2n-1)√(1+Q^2π ^2/β ^2S^2) 
    ] (4S)^n/n-1

In order to check the validity of the first law of thermodynamics for our
solutions, we consider entropy S and electric charge Q as a complete set
of extensive quantities for mass. Also, we define temperature T and
electric potential U as their conjugate intensive quantities, respectively:

    T=( ∂ M(S,Q)/∂ S) _Q    and    U=( ∂ M(S,Q)/∂ Q) _S.

Our numerical calculations show that, above quantities satisfy the first law
of black hole thermodynamics: 

    dM=TdS+UdQ.

In following section, we will investigate the thermal stability of obtained
solutions.



§ THERMAL STABILITY 




This section is devoted to study the thermal stability of the obtained
solutions. One can consider the stability of a black hole in both canonical
and grand canonical ensembles. In canonical ensemble, positivity of the heat
capacity together with the positive values for the temperature guaranties
thermal stability of the solutions. Heat capacity is defined as 

    C_Q=T( ∂ S/∂ T) _Q=( 
    ∂ M/∂ S) _Q( ∂ ^2M/∂
    S^2) _Q^-1.

In case of unstable black hole, a phase transition occur to get stable
state. We desire to investigate the phase transition point by finding the
roots and divergencies of C_Q. The roots of heat capacity which are
equivalent to zero values of temperature, determine a thermal transition
between un-physical (T<0) and physical (T>0) state of the black hole.
Beside, we recall that based on the thermal physics, when the heat capacity
diverges, a second order phase transition will occur.

Since it is not easy to determine the roots and divergencies of the heat
capacity analytically, we study this by using the figures. The behaviors of C_Q in terms of r_+ for different sets of parameters have been
depicted in Figs. <ref>-<ref>. It is obvious that the stability
depends on the metric parameters. By choosing a proper set of metric
parameters, Fig. <ref> specifies that by increasing the electric charge 
Q (massive parameter m), the transition between un-physical and physical
states will occur for larger black holes.

Fig. <ref> illustrates a second order phase transition of our solutions
for some specific values of the model parameters. The transition is between
two physical states since the temperature is positive. In Fig. <ref>,
we observe that for smaller black holes C_Q is negative, so they are not
stable and transition allows us to have larger stable black holes. In
addition, as we decrease the nonlinearity of electrodynamics (larger β
), the transition happens for smaller black holes (Fig. <ref>). As Fig. 
<ref> shows, for higher dimensions phase transition occurs for larger
horizon radii.

Changing massive model parameters causes rich phase transitions (Fig. <ref>). The solid curve in Fig. <ref> shows two transitions. One from
stable small black hole to unstable middle black hole and another from
unstable middle black hole to stable large black hole. These transitions
represent that depending on the metric parameters, the middle black holes
may not be allowed. This fact is clear from dash-double-dotted curve for
which the transition is just from small to large black holes. 



Now, we turn to study the stability in the grand canonical ensemble.
Basically, the local stability can be carried out by finding the determinant
of the Hessian matrix 𝐇 with respect to its extensive variables, X_i's. In the canonical ensemble, Q is a fixed parameter so the
positivity of heat capacity (C_Q) is sufficient to ensure the local
thermal stability. However, in grand canonical ensemble, the number of
extensive parameters depends on the theory, and we have to check the sign of 
𝐇 to search for thermal stability. In our case, the mass is a
function of extensive variables entropy S and charge Q. So, Hessian
matrix defines as


    𝐇_X_i,X_j^M=∂ ^2M/∂ X_i∂
    X_J=[ 
    [  ∂ ^2M/∂ S^2 ∂ ^2M/∂
    S∂ Q; ∂ ^2M/∂ S∂ Q  ∂ ^2M/∂
    Q^2 ]
    ].

It is not easy to investigate the thermal stability analytically, so we turn
to the figures again. The behavior of 𝐇 in terms of r_+ for
different sets of model parameters has been plotted in Fig. <ref> and 
<ref>. By fixing other metric parameters, Fig.<ref> shows a
minimum value for horizon radius that for values under which black holes are
unstable. This minimum value grows with increasing the mass parameter of
massive gravity. The same behavior could be seen as well in terms of
dimension n (Fig. <ref>). We emphasize that in this range the
temperature is positive.

Changing nonlinear parameter β, may create two roots for 𝐇
(Fig. <ref>). But since the temperature is negative for smaller root,
the only important limit for allowed horizon radius value depends on the
larger root. As one can see from Fig. <ref>, the horizon of stable
black hole should be larger than larger root. It is notable to mention that
by increasing β this limit is removed. We observe from Fig.<ref> that increasing electric charge Q causes two roots for 𝐇. We do not care to the smaller root since it places in the negative range
of temperature. Furthermore, the larger root grows for smaller electric
charges, so the stable black holes have larger horizons.



§ CRITICAL BEHAVIOR AND REENTRANT PHASE TRANSITION 








In this section, we first explore critical behavior of (3+1)-dimensional
logarithmic charged massive black hole solutions with spherical event
horizon (k=1) by investigating the behavior of the specific heat at
constant charge, C_Q. Then, by studying Gibbs free energy G, we seek
for different phase transitions of systems.

The specific heat at constant charge is defined as

    C_Q=T( dS/dT) _Q.

The sign of C_Q indicates the local thermodynamic stability/instability
of the system. The positive (negative) sign of this quantity shows the local
thermodynamic stability (instability). Note that Eq. (<ref>) is also
calculated at fixed cosmological constant and parameter of nonlinearity β. In Fig. <ref>, we plot the behavior of temperature with
respect to event horizon radius for different values of the charge. The sign
of slope in T-r_+ curves in allowed regions with positive temperature
determines the sign of C_Q. Expanding the Hawking temperature (<ref>
) for small r_+, we receive

    T=4β/r_+( Q_m-Q) +c_0c_1m^2/4π+
    r_+( 3+4l^2β ^2[ 1+ln( 2π Q/β) 
    ] ) /4π l^2-2β ^2r_+ln( r_+) 
    /π+O( r_+^3) ,

where Q_m=( 1+c_0^2c_2m^2) /( 16πβ) is the `marginal charge' and we use Λ =-3/l^2. As one
can see, the behavior of the temperature for small r_+ depends on the
charge of the black hole. For Q>Q_m, black hole has zero temperature for
a specific value of r_+ and for values greater than it temperature is
positive. It means that these black holes are Reissner-Nordstrom (RN) type.
For Q<Q_m, the black hole does not exist in low temperature region and
temperature diverges as r_+ goes to zero (similar to Schwarzschild
solution). So, these black holes are Schwarzschild (S) type. Note that for Q=Q_m, the black hole has a finite positive temperature as r_+ tends
to zero. Therefore, the latter is neither S-type, nor RN-type and is
marginal. This is the effect of massive term c_0c_1m^2/4π in (<ref>). For massless case (m=0), T→ 0 as r_+→
0 for Q=Q_m and the black hole is RN-type <cit.>. The large 
r_+ limit of the temperature is 3r_+/4π l^2 which is clearly
independent of the charge and increases linearly as r_+ increases (see
Fig. <ref>).

Now, we intend to obtain the critical point (where continuous phase
transition occurs). In order to characterize the critical point, we have

    ∂ T/∂ S|_Q_c=0,      
    ∂ ^2T/∂ S^2|_Q_c=0,

where Q_c is the critical charge. Using (<ref>), one could
determine Q_c and r_+c and specify the critical point (Q_c,T_c,r_+c). In Fig. <ref>, we plot the behavior of Q_c
, T_c and r_+c with respect to β for a set of model
parameters by solving the equations of (<ref>) numerically. The
vertical dashed lines in Fig. <ref> show the range in which the
critical behavior takes place for S-type black holes. For β values
larger than one specified by the rightmost vertical dashed line, the
critical point occurs for RN-type black holes.

Now, we turn to study the behavior of Gibbs free energy for our black holes
to find out the phase structure of the system. It determines the globally
stable state at equilibrium process. The lowest Gibbs free energy specifies
the global stable state. One could obtain the Gibbs free energy per unit
volume V_2 as <cit.>

    G=M-TS,

We focus on RN-type (3+1)-dimensional logarithmic charged massive black
holes. The Gibbs free energy as a function of temperature for different
values of charge is illustrated in Fig. <ref>. In addition, the
phase diagram is depicted in the Q–T diagram in Figs. <ref>.
According to the value of charge Q, we have different phase transitions.
For Q=Q_c, the system is at critical point. This point is highlighted by
a black spot in Fig. <ref>. For charge values greater than Q_c,
the system is locally stable (C_Q>0) and Gibbs free energy is single
valued. Therefore, no phase transition happens. For charge values less than
critical value Q_c, there are thermally instable phases in the system (C_Q<0). So, the system experiences different phase transitions depending
on the value of charge Q. A first order phase transition between small and
large black holes occurs for Q_z<Q<Q_c. This is accompanied by a
discontinuity in the slope of Gibbs free energy at transition point. This
phase transition is marked by the gold curve in Fig. <ref>. For Q
in the range of Q_t<Q<Q_z, an interesting phenomenon happens. In this
range, in addition to a first order phase transition which sepatates small
and large black holes, a zeroth order phase transition between small and
intermediate black holes occurs. This can be seen as a finite jump in Gibbs
free energy in Fig. <ref>. The latter initiates from Q=Q_z and
terminates at Q=Q_t. z and t points are specified in Fig. <ref>. In this figure, this zeroth order phase transition is identified
by purple curve. This intermediate (large)/small/large transition shows a
reentrant phase transition <cit.>. This is the first observation
of reentrant phase transition for black hole solutions in massive gravity
without extending thermodynamics phase space. In <cit.>, the occurance
of reentrant phase transition in extended thermodynamics phase space has
been reported for black holes in massive gravity. As Fig. <ref>
exhibits, no phase transition happens for Q<Q_t and only stable large
black holes exists. Also, no black hole region in Fig. <ref>
corresponds to no black hole solution exists.



§ SUMMARY AND CLOSING REMARKS


The importance of massive gravity model is at least twofold. On the one
hand, because of the recent developments in gravitational waves observatory
indicating the speed of gravitational wave deviates from the speed of light
which implies the existence of massive graviton, we need an extension of
general relativity with a massive graviton as intermediate particle similar
to W and Z in electroweak interaction. On the other hand, from
holographic viewpoint, giving mass to gravity could model momentum
dissipation in dual system on boundary. In this work, we have studied black
hole solutions in massive gravity model in the presence of LNE Lagrangian.
As well as BI model, in the logarithmic one, there is no divergency in the
electric field of a point particle at the origin.

We first constructed the field equations by varying the action. Then, we
obtained the exact solution of the metric function and studied the behavior
of the obtained solutions. We observed that depending on the model
parameters, it is possible to have two horizons, an extreme horizon, or
naked singularity. Moreover, one of the effects of massive gravity is the
existence of multi horizons. This phenomenon is due to anti-evaporation
which is a quantum effect. Next, we calculated the conserved and
thermodynamic quantities. Calculation of these quantities is necessary for
checking the satisfaction of the first law of black hole thermodynamics. It
was observed that in the presence of nonlinear electrodynamics and within
the context of massive gravity model, the first law of thermodynamic is
still valid.

We studied the thermal stability in both canonical and grand canonical
ensembles. We observed that considering different values for the model
parameters leads to a variety of phase structures for our black holes. In
canonical ensemble, the heat capacity can experience both positive and
negative values. The position of its root(s) and divergency(ies) affected by
metric parameters. In this ensemble, our black holes encounter two types of
phase transitions; the second order phase transition in the position of
divergences of heat capacity with positive temperature, and transitions
between un-physical (T<0) and physical (T>0) states of solutions.
Increasing the values of electric charge (Q) or massive parameter (m)
causes the transition between un-physical and physical black holes in larger
sizes (Fig. <ref>). It means that there is a lower bound on black hole
radius increased as Q and m are enhanced. Also, by decreasing the
nonlinearity parameter (β), second order phase transition happens for
smaller black holes (Fig. <ref>). In the grand canonical ensemble, the
determinant of the Hessian matrix 𝐇 had a root with positive
temperature, which its value depends on the metric parameters. We found out
that some horizon radii are not allowed for our solution since 𝐇
is negative for them. For instance, when other parameters are fixed, the
stability is observed for larger black hole when massive parameter m (or
dimension n) increases (Fig. <ref>). Also, for smaller values of
electric charge, the stability is guaranteed for larger black holes (Fig. 
<ref>).

Finally, we investigated the critical behavior of (3+1)-dimensional black
holes with spherical event horizon. We showed that there is a marginal
charge Q_m that for charge values greater than it black holes are
Reissner-Nordstrom-type, whereas for charge values less than it black holes
are Schwarzschild-type. For massless gravity, the black hole solutions with Q=Q_m are RN-type. We observed that in the context of massive gravity the
latter are neither S-type nor RN-type and have a marginal behavior due to
the effects of massive term. Next, we focused on RN-type black holes and
study the behavior of Gibbs free energy to disclose the phase structure of
solutions. We found out that there are different phase transitions according
to the value of charge Q. There are some ranges for which there is no
phase transition. However, for another range we just have a first order
phase transition between small and large black holes accompined by a
discontinuity in the slope of Gibbs free energy at transition point.
Interestingly enough, there is a range of charge values for which black
holes experience a zeroth order phase transition where a finite jump in
Gibbs free energy value occurs. This causes a reentrant phase transition
i.e. a large/small/large phase transition. This is the first report of
reentrant phase transition for black holes in massive gravity without
needing to extend thermodynamics phase space.

The study of this work can be extended from differnt aspects. For instance,
in this context, one could explore the behavior of holographic conductivity.
Furthermore, the effects of massive gravity with logarithmic electrodynamics
on holographic superconductors could be investigated. Here, we examined the
critical behavior of (3+1)-dimensional black holes numerically. This may
be studied analytically. Also, one could explore this behavior for higher
dimensional black holes. Moreover, the effects of more model parameters on
phase behavior of such systems could be taken under study. Some of these
pointed out issues are now under investigation.


AS and SMZ thank the Research Council of Shiraz University. MKZ would like
to thank Shahid Chamran University of Ahvaz for supporting this work. We
gratefully acknowledge useful discussions with A. Dehyadegari.


99
ref1 S. N. Gupta, Phys. Rev. 96, 1683 (1954);S. Weinberg, Phys. Rev. B 988, 138 (1965);S. Deser, Gen. Relativ. Gravit. 1, 9 (1970).

Abbot1 B. P. Abbott et al. [LIGO Scientific and Virgo
Collaborations], Phys. Rev. Lett. 119, 161101 (2017).

Abbot2 B. P. Abbott et al.[LIGO Scientific and Virgo
Collaborations], Phys. Rev. Lett. 116, 061102 (2016).

Ales A. Casalino, M. Rinaldi, L. Sebastiani, S. Vagnozzi, Phys.
Dark Univ. 22, 108 (2018).

ref8 D. Vegh, arXiv:1301.0537 [hep-th].

ref3 M. Fierz and W. Paul, Proc. R. Soc. A 173, 211(1939).

ref4 H. van Dam and M. J. G. Veltman, Nucl. Phys. B 22,
397 (1970).

ref5 D. G. Boulware and S. Deser, Phys.Rev. D 6, 3368
(1972).

ref16 A. I. Vainshtein, Phys. Lett. B 39, 393 (1972).

ref6 C. de Rham, G. Gabadadze and A. J. Tolley, Phys. Rev. Lett.,
106, 231101 (2011).

ref13 A. H. Chamseddine and V. Mukhanov, JHEP 1303, 092
(2013);S. Deser and A. Waldron, Phys. Rev. Lett. 110 ,111101 (2013).

ref11 A. Salam and J.A. Strathdee, Phys. Rev. D 16, 2668
(1977);Z. Berezhiani, D. Comelli, F. Nesti and L. Pilo, JHEP 07, 130
(2008);T.M. Nieuwenhuizen, Phys. Rev. D 84, 024038 (2011);L. Berezhiani, G. Chkareuli, C. de Rham, G. Gabadadze and A. J. Tolley,
Phys. Rev. D 85, 044024 (2012);R. Brito and V. Cardoso, P. Pani, Phys. Rev. D 88, 064006 (2013);
E. Babichev and A. Fabbri, JHEP 1407, 016 (2014);A. Adams, D.A. Roberts and O. Saremi, Phys. Rev. D 91(4), 046003
(2015);J. Xu, L.M. Cao, Y.P. Hu, Phys. Rev. D 91, 124033 (2015).

ref7 Th. M. Nieuwenhuizen, Phys. Rev. D 84, 024038 (2011);
K. Koyama, G. Niz, and G. Tasinato, Phys. Rev. D 84, 064033 (2011);
L. Berezhiani, G. Chkareuli, C. de Rham, G. Gabadadze, and A. J. Tolley,
Phys. Rev. D 85, 044024 (2012).

ref9 M. Blake and D. Tong, Phys. Rev. D 88, 10, 106004
(2013);R. A. Davison, Phys. Rev. D 88, 086003 (2013);R. A. Davison, K. Schalm and J. Zaanen, Phys. Rev. B 89, 245116
(2014).

ref10 R. G. Cai, Y.P. Hu, Q. Y. Pan and Y. L. Zhang, Phys. Rev. D 
91, 024032 (2015).

ref12 A. Ghodsi and D. Mahdavian, JHEP, 1206, 131 (2012).

ref17 A. Gruzinov and M. Mirbabayi, arXiv:1106.2551 [hep-th].

ref18 T. Katsuragawa, S. Nojiri, S.D. Odintsov and M. Yamazaki,
Phys. Rev. D 93, 124013 (2016).

cosmicAcc G. Dvali, G. Gabadadze and M. Shifman, Phys. Rev. D 
67, 044020 (2003);G. Dvali, S. Hofmann and J. Khoury, Phys. Rev. D 76, 084006 (2007);
B. Eslam Panah and S. H. Hendi, EPL 125, 60006 (2019).

ref14 C. de Rham, Living Rev. Relativity 17, 7 (2014).

ref15 K. Hinterbichler, Proceedings of the 51st Rencontres de
Moriond, ARISF, ISBN: 979-10-968-7901-4, (2016).

ref19 F. Capela and P. G. Tinyakov, JHEP 1104, 042 (2011).

ref20 M. Born and L. Infeld, Proc. R. Soc. Lond. 144, 425
(1934).

ref21 W. Heisenberg and H. Euler, Z. Phys. 98, 714 (1936).

ref22 J. Plebanski, [A Study], C.I.E.A. del I.P.N., Mexico City,
(1966).

ref23 N. Seiberg and E. Witten, JHEP 09, 032 (1999).

ref26 E. Fradkin and A. Tseytlin, Phys. Lett. B 163, 123
(1985);E. Bergshoeff, E. Sezgin, C. Pope and P. Townsend, Phys. Lett. B 188
, 70 (1987);R. Metsaev, M. Rahmanov and A. Tseytlin, Phys. Lett. B 193, 207
(1987);C. Callan, C. Lovelace, C. Nappi and S. Yost, Nucl. Phys. B 308,
221 (1988).

ref24 M. Halilsoy, O. Gurtug and S. H. Mazharimousavi, Astropart.
Phys. 68, 1 (2015).

ref25 J. Bardeen, presented at GR5, Tiflis, U.S.S.R., and
published in the conference proceedings in the U.S.S.R. (1968).

ref27 H. Soleng, Phys. Rev. D 52, 6178 (1995).

ref28 S. H. Hendi, JHEP 03, 065 (2012);S. H. Hendi, Ann. Phys. 333, 282 (2013).

ref30 A. Sheykhi, F. Naeimipour and S. M. Zebarjad, Phys. Rev. D 
91, 124057 (2015).

ref31 S. Hajkhalili and A. Sheykhi, Inter. J. Mod Phys. D 
27, 1850075 (2018).

ref32 A.Dehyadegari, A.Sheykhi and M. Kord Zangeneh, Phys. Lett. B 
758, 226 (2016).

ref33 A. Sheykhi and Z. Mahmoudi, Gen. Relativ. Gravit. 47
, 90 (2015).

ref34 A. Sheykhi, F. Naeimipour and S. M. Zebarjad, Gen. Relativ.
Gravit. 48, 33 (2016).

ref35 S. H. Hendi, B. Eslam Panah and S. Panahiyan, Phys. Rev. D 
91, 084031 (2015).

ref38 S. H. Hendi, JHEP 1203, 065 (2012);S. H. Hendi, Ann. Phys. 346, 42 (2014);A. Sheykhi and F. Shamsi, Int J Theor Phys 56, 916 (2017);Z. Dayyani, A. Sheykhi, M. H. Dehghani and S. Hajkhalili, Eur. Phys. J. C 
78, 152 (2018);A. Sheykhi, D. Hashemi Asl and A. Dehyadegari, Phys. Lett. B 781,
139 (2018);A. Dehyadegari, A. Sheykhi and M. Kord Zangeneh, Phys. Lett. B 758,
226 (2016).

maldacena J. M. Maldacena, Adv. Theor. Math. Phys. 2, 231
(1998).

refmetric R. G. Cai, Y. P. Hu, Q. Y. Pan and Y. L. Zhang, Phys.
Rev. D 91, 024032 (2015);J. Xu, L. M. Cao and Y. P. Hu, Phys. Rev. D 91, 124033 (2015);
S. H. Hendi, B. Eslam Panah and S. Panahiyan, JHEP 1511, 157 (2015).

ref39 J. D. Beckenstein, Phys. Rev. D 7, 2333 (1973);
S. W. Hawking and C. J. Hunter, Phys. Rev. D 59, 044025 (1999).

ref40 T. Katsuragawa, Universe 1, 158 (2015);T. Katsuragawa and S. Nojiri, Phys. Rev. D 91, 084001 (2015).

ref41 R. G. Cai, Y. P. Hu, Q. Y. Pan and Y. L. Zhang, Phys. Rev. D 
91, 024032 (2015);J. Xu, L. M. Cao and Y. P. Hu, Phys. Rev. D 91, 124033 (2015).

1711.01151 A. Dehyadegari and A. Sheykhi, Phys. Rev. D 98
, 024011 (2018).

RPT1 S. Gunasekaran, D. Kubiznak and R. B. Mann, JHEP 1211,
 110 (2012).

RPT2 N. Altamirano, D. Kubiznak and R. B. Mann, Phys. Rev. D 
88, 101502 (2013).

RPT3 D. C. Zou, R. Yue and M. Zhang, Eur. Phys. J. C 77,
256 (2017);M. Zhang, D. C. Zou and R. Yue, Advances in High Energy Physics 2017
, Article ID 3819246 (2017).


