# Charged $AdS$ Black Holes in $4D$ Einstein–Gauss–Bonnet Massive Gravity

Prosenjit Paul<EMAIL_ADDRESS>Indian Institute Of Engineering
Science and Technology (IIEST), Shibpur, WB 711103, India  Sudhaker
Upadhyay111Corresponding author222Visiting Associate, IUCAA Pune,
Maharashtra-411007, India<EMAIL_ADDRESS>Department of Physics,
K. L. S. College, Magadh University, Nawada, Bihar 805110, India  School of
Physics, Damghan University, Damghan, 3671641167, Iran  Dharm Veer
Singh333Visiting Associate, IUCAA Pune, Maharashtra-411007, India.
<EMAIL_ADDRESS>Department of Physics, Institute of Applied Science and
Humanities, GLA University, Mathura, 281406 India

###### Abstract

We investigate Einstein–Gauss–Bonnet–Maxwell massive gravity in $4D$ AdS
background and find an exact black hole solution. The horizon structure of the
black holes studied. Treating the cosmological constant as pressure and Gauss-
Bonnet coupling parameters, and massive gravity parameters as variables, we
drive the first law of black hole thermodynamics. To study the global
stability of the black holes we compute the Gibbs free energy. The local
stability of the black hole is also studied through specific heat. We analyze
the effects of graviton mass and Gauss-Bonnet coupling parameters on the phase
transition of the black holes. Finally, the effects of graviton mass and
massive gravity parameters on the Joule-Thomson expansion of the black hole
are studied.

##  I Introduction

General Relativity (GR) is a theory of gravitation that help us understand
gravitational waves, gravitational lensing, an effect of gravity on time known
as gravitational time dilation and black holes. Although GR is not a complete
theory of quantum gravity, the simplest theory of gravity describes various
astrophysical phenomena. To complete the GR, people try to modify it in
various ways, for instance, by adding the higher order term to the Einstein-
Hilbert action of GR but still, a complete theory is missing. Some other
examples of higher order gravity theory are scalar-tensor theories [1, 2, 3,
4, 5, 6], Lovelock gravity [7, 8, 9], regular black holes [10, 11, 12, 13] and
brane world cosmology [14, 15, 16].

Lovelock theories [7, 8] is a special theory of higher-order gravity in $4D$
spacetime that preserves diffeomorphism invariance, metricity and second-order
equations of motion. From Lovelock’s theories of gravity, Gauss-Bonnet gravity
can be obtained in higher dimensions [17]. The Gauss-Bonnet term does not
contribute to the dynamics of the theory in four dimensions but rather
contributes to the dynamics when the dimensions of spacetime are greater than
four. In recent days, Glavan and Lin [18] found the solution to the Einstein-
Gauss–Bonnet field equation in four dimensions by rescaling Gauss-Bonnet
coupling parameter $\alpha$ by $\alpha/D-4$. However, the charged AdS solution
of Einstein-Gauss–Bonnet theory was found in Ref. [19]. For a complete
discussion on $4D$ Gauss-Bonnet gravity, see Refs. [20]. Some other static
spherically symmetric black hole solutions and their thermodynamics, phase
transition in $4D$ or higher dimensions Einstein-Gauss–Bonnet gravity studied
in Refs. [21, 23, 22, 24, 25, 26, 27, 28]. Einstein-Gauss–Bonnet black hole
solution in nonlinear electrodynamics is studied in Refs. [29, 30, 32, 35, 33,
36, 37, 34, 31].

Another way of modifying the GR is by adding a mass to the graviton. According
to the GR, the graviton is a massless spin-2 particle. But one may ask if that
is a self- a consistent theory of massive gravity possible or not. In fact,
people have tried to answer this question by modifying the Einstein-Hilbert
action that describes massive graviton. The recent observation of
gravitational waves by LIGO put a maximum limit on the graviton mass, $m\leq
1.2\times 10^{-22}$eV [38]. A theory of massive gravity was first constructed
by Fierz and Pauli in 1939 [39, 40]. In the curved background, this theory
encounters ghost instabilities [41]. A new nonlinear massive gravity theory
proposed by de Rham, Gabadadze and Tolley (dRGT) [42, 43] avoids ghost
problems. The charged black holes in Gauss–Bonnet massive gravity studied
[44]. Some other spherically symmetric black holes in massive gravity and
their thermodynamics have also been studied [45, 48, 49, 46, 47, 50].

In our theory, we consider the anti-de-Sitter (AdS) background, i.e. we add a
negative cosmological constant. The negative cosmological constant is a
crucial ingredient in the AdS/CFT correspondence, a duality between a theory
of quantum gravity in (AdS) space and a conformal field theory (CFT) in one
lower dimension. The AdS/CFT correspondence allows for the study of strongly
coupled field theories using classical gravity, which is a useful tool for
investigating non-perturbative phenomena that cannot be understood through
standard perturbative methods. One significant consequence of the negative
cosmological constant is that it leads to the presence of a holographic screen
at the AdS boundary, which encodes the bulk geometry’s information. This
holographic principle means that the number of degrees of freedom in the AdS
space is proportional to the area of the holographic screen, rather than the
volume, as in ordinary theories. The AdS/CFT correspondence, therefore,
implies that the degrees of freedom in the AdS space are equivalent to those
of the boundary CFT. The negative cosmological constant also plays a critical
role in the AdS black hole physics. Black holes in AdS space can have a
negative specific heat, which is impossible in flat space. This phenomenon is
related to the AdS space’s boundary conditions, which force the black hole to
lose energy and mass through the AdS boundary, leading to a reduction in
temperature. This behavior is known as Hawking-Page phase transition [51],
where the black hole is in thermal equilibrium with a thermal AdS space. The
AdS/CFT correspondence allows for the study of the thermodynamic behaviour of
black holes using the corresponding CFT, providing insights into the nature of
black hole thermodynamics.

Recently, researchers have considered the cosmological constant as a variable
parameter and linked it to the thermodynamic pressure, which is conjugate to
the thermodynamic volume [52, 53, 54, 55, 56]. This approach has resulted in
an extended phase space, where the black hole mass is regarded as the
enthalpy, instead of the internal energy [53]. Many studies have explored the
thermodynamics and phase transitions of black holes in this extended phase
space, revealing new phenomena such as P-V criticality in various black holes
spacetime [57, 58].

In AdS spacetime, the black hole mass is naturally treated as the enthalpy,
leading to the consideration of Joule-Thomson expansion for the black hole.
This expansion investigates isenthalpic curves, which are constant mass
curves. Previous investigations of the Joule-Thomson expansion for charged AdS
[59] and Kerr AdS [60] black holes have been conducted within the framework of
Einstein gravity. However, extended theories of Einstein gravity introduce new
physical degrees of freedom, raising questions about their role and physical
impact on the Joule-Thomson expansion. In this paper, we investigate how the
presence of massive gravity modifies the Joule-Thomson expansion of the
charged AdS black hole in Gauss–Bonnet gravity, inspired by recent progress in
understanding massive gravity. The approach taken here is relevant not only
for charged AdS black holes in Eiantein-Gauss-Bonnet massive gravity but also
for those in other alternative theories of gravity where additional
gravitational modes emerge. For instance, the Joule-Thomson expansion may be
examined for charged AdS black holes in teleparallel $f(T)$ gravity [61, 62],
with $T$ representing torsion, or in $f(R)$ gravity [63] with a nonlinear
electrodynamics field.

Since the charged AdS Einstein-Gauss–Bonnet theory in massive gravity is not
studied yet. Therefore, in this paper, we investigate charged Einstein-
Gauss–Bonnet massive gravity and find an exact solution in $4D$ AdS space. we
also discuss the horizon structure of charged AdS black hole in $4D$
Einstein–Gauss–Bonnet massive gravity. Moreover, we discuss the thermal
properties of this black hole. To be more precise, we compute the entropy and
temperature that satisfy the first law of black hole thermodynamics. In AdS
space, the mass of the black hole is treated as enthalpy. Furthermore, we
analyzed the stability and Van der walls like phase transition of the black
holes. Here, we study the Joule-Thomson expansion of charged AdS black hole in
4D Einstein–Gauss–Bonnet massive gravity, and the constant mass curves are
known as isenthalpic curves. Finally, we investigate the effects of graviton
mass and the massive gravity parameters on the Joule-Thomson expansion of
charged AdS black holes in $4D$ Einstein–Gauss–Bonnet massive gravity.

The paper is organized as follows. In section II, we discuss the action
describing the Einstein–Gauss– Bonnet–Maxwell massive gravity in $4D$ AdS
space and their field equations. Here, we find the exact black hole solution.
The effects of graviton mass on the horizon structure of the black hole are
also depicted. In section III, we study the first law of black hole
thermodynamics and the effects of graviton mass on Hawking temperature. To
investigate the global stability of the black holes, we compute the Gibbs free
energy. Next, in section V, the effects of graviton mass on the local
stability of the black hole are studied. The Van der walls-like phase
transition of the black hole is analyzed in section VII. We numerically
investigate the effects of the graviton mass, the charge of the black hole and
Gauss–Bonnet coupling parameter on the critical parameters (namely, critical
volume, critical pressure and critical volume) of the black hole. The effects
of critical parameters on the phase transition of the black hole are also
studied. we investigate the Joule-Thomson expansion of the black hole in
section VIII. Here, we analyze the effects of graviton mass and massive
gravity parameters on the constant mass curve and inverse curve. Finally, we
compute Joule-Thomson thermodynamic coefficient as a function of the black
hole horizon radius.

##  II Einstein-Gauss–Bonnet Massive Gravity in 4D

The action for Einstein–Maxwell–Gauss–Bonnet massive gravity with a negative
cosmological constant in $D$ dimensions is given by

$S=\frac{1}{16\pi}\int
d^{D}x\sqrt{-g}\Biggr{[}R-2\Lambda+\alpha\mathcal{G}-F_{\mu\nu}F^{\mu\nu}+m^{2}\sum_{i}c_{i}\mathcal{U}_{i}(g,h)\Biggr{]},$
(1)

where $g$ is determinant of the metric $g_{\mu\nu}$, $R$ is Ricci scalar,
$\alpha$ is Gauss–Bonnet coupling parameter,
$\mathcal{G}=R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}-4R_{\mu\nu}R^{\mu\nu}+R^{2}$
is the Gauss–Bonnet term, $R_{\mu\nu\rho\sigma}$ is Riemann tensor,
$R_{\mu\nu}$ is Ricci tensor and
$F_{\mu\nu}=\partial_{\mu}{A_{\nu}}-\partial_{\nu}{A_{\mu}}$ is Maxwell
tensor. Apart from that $m$ is a parameter related to graviton mass,
$h_{\alpha\nu}$ is a fixed symmetric tensor and usually is called the
reference metric, $c_{i}(i=1,2,3,4)$ are constant444In order to have a self-
consistent massive gravity theory, the coupling parameters $c_{i}$ might be
required to be negative if the squared mass of the graviton is positive.
However, in the AdS spacetime, the coupling parameters $c_{i}$ can still take
the positive values. This is because the fluctuations of the fields with the
negative squared masses in the AdS spacetime could still be stable if their
squared masses obey the corresponding Breitenlohner–Freedman bounds. [64] and
$\mathcal{U}_{i}(g,h)$ is symmetric polynomials of eigenvalues of matrix
$\mathcal{K}_{\nu}^{\mu}=\sqrt{g^{\mu\alpha}h_{\alpha\nu}}$, given by

$\begin{split}\mathcal{U}_{1}&=\bigr{[}\mathcal{K}\bigr{]},\\\
\mathcal{U}_{2}&=\bigr{[}\mathcal{K}\bigr{]}^{2}-\bigr{[}\mathcal{K}^{2}\bigr{]},\\\
\mathcal{U}_{3}&=\bigr{[}\mathcal{K}\bigr{]}^{3}-3\bigr{[}\mathcal{K}\bigr{]}\bigr{[}\mathcal{K}^{2}\bigr{]}+2\bigr{[}\mathcal{K}^{3}\bigr{]},\\\
\mathcal{U}_{4}&=\bigr{[}\mathcal{K}\bigr{]}^{4}-6\bigr{[}\mathcal{K}^{2}\bigr{]}\bigr{[}\mathcal{K}\bigr{]}^{2}+8\bigr{[}\mathcal{K}^{3}\bigr{]}\bigr{[}\mathcal{K}\bigr{]}+3\bigr{[}\mathcal{K}^{2}\bigr{]}^{2}-6\bigr{[}\mathcal{K}^{4}\bigr{]},\end{split}$
(2)

where parentheses $[...]$ represents trace of the matrix
$\mathcal{K}_{\nu}^{\mu}$. In $D=4$ dimensions the Gauss–Bonnet term does not
contribute to the dynamics, so we rescale the Gauss–Bonnet coupling parameter
$\alpha\to\alpha/(D-4)$ [18], therefore, the action takes the following form:

$S=\frac{1}{16\pi}\int
d^{D}x\sqrt{-g}\left[R-2\Lambda+\frac{\alpha}{D-4}\mathcal{G}-F_{\mu\nu}F^{\mu\nu}+m^{2}\sum_{i}c_{i}\mathcal{U}_{i}(g,h)\right].$
(3)

Now, we consider a static and spherically symmetric solution of the form

$ds^{2}=-e^{2A(r)}dt^{2}+e^{2B(r)}dr^{2}+r^{2}d{\Omega}_{D-2}^{2},$ (4)

and following the Ref. [65] we take the reference metric as

$h_{\mu\nu}=diag\bigl{(}0,0,c^{2},c^{2}\sin^{2}\theta\bigl{)},$ (5)

where $c$ is a dimensionless positive constant. The reference metric
$h_{\mu\nu}$ is a rank two symmetric tensor. Physically, $h_{\mu\nu}$
corresponds to the background metric around which fluctuations take the
Fierz–Pauli form. By using equations (2) and (5), we obtain

$\begin{split}\mathcal{U}_{1}&=\frac{(D-2)c}{r},\\\
\mathcal{U}_{2}&=\frac{(D-2)(D-3)c^{2}}{r^{2}},\\\
\mathcal{U}_{3}&=\frac{(D-2)(D-3)(D-4)c^{3}}{r^{3}},\\\
\mathcal{U}_{4}&=\frac{(D-2)(D-3)(D-4)(D-4)c^{4}}{r^{4}}.\end{split}$ (6)

Substituting the metric and the electrostatic potential in action (3), the
first integral exists [19]

$\phi^{\prime}(r)=-\frac{Q}{r^{D-2}}e^{A+B},$ (7)

and taking the limit $D\to 4$ and using the relation $\Lambda=-3/l^{2}$ we
obtain

$S=\frac{\Sigma_{2}}{16\pi}\int
dtdr2e^{A+B}\biggr{[}r^{3}\psi\Bigl{(}1+\alpha\psi\Bigl{)}+\frac{r^{3}}{l^{2}}+\frac{Q^{2}}{r}+m^{2}\Bigl{\\{}\frac{c_{1}cr^{2}}{2}+{c_{2}c^{2}r}\Bigl{\\}}\biggr{]}^{\prime},$
(8)

where prime denotes differentiation with respect to r,
$\Sigma_{2}=\frac{2{\pi}^{\frac{3}{2}}}{\Gamma\Bigl{(}{1+\frac{1}{2}}\Bigl{)}}$
and $\psi=r^{-2}\Bigl{(}1-e^{-2B}\Bigl{)}$ with

$e^{A+B}=1.$ (9)

If we choose $m=0$ or $c=0$ then equation (8) reduced to action in massless
gravity. Now, using the action (8), we obtain solution as

$\psi\Bigl{(}1+\alpha\psi\Bigl{)}+\frac{1}{l^{2}}+\frac{Q^{2}}{r^{4}}+\frac{m^{2}}{r^{3}}\Bigl{\\{}\frac{c_{1}cr^{2}}{2}+{c_{2}c^{2}r}\Bigl{\\}}-\frac{8\pi
M}{\Sigma_{2}r^{3}}=0,$ (10)

where $M$ is the integration constant related to the mass of the black hole.
Therefore, the exact solution is

$e^{2A}=e^{-2B}=1+\frac{r^{2}}{2\alpha}\Biggr{[}1\pm\sqrt{1+4\alpha\biggl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}}{2r^{2}}\Bigl{(}cc_{1}r+2c^{2}c_{2}\Bigl{)}\biggl{\\}}}\Biggr{]}.$
(11)

The negative branch corresponds to the $4D$ charged AdS EGB massive black
hole, whereas the +ve branch does not lead to a physically meaningful solution
because the positive sign in the mass term indicates graviton instabilities,
so we only take the negative branch of equation (11). For the chargeless
limit, our solution reduces to a black hole solution as obtained in Ref. [50]

$e^{2A}=e^{-2B}=1+\frac{r^{2}}{2\alpha}\Biggr{[}1-\sqrt{1+4\alpha\biggl{\\{}\frac{2M}{r^{3}}-\frac{1}{l^{2}}-\frac{m^{2}}{2r^{2}}\Bigl{(}cc_{1}r+2c^{2}c_{2}\Bigl{)}\biggl{\\}}}\Biggr{]}.$
(12)

In the limit $\alpha\to 0$, equation (11) reduces to the charged AdS black
hole in massive gravity [66]

$e^{2A}=e^{-2B}=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}+\frac{r^{2}}{l^{2}}+\frac{m^{2}}{2}\bigl{(}cc_{1}r+2c^{2}c_{2}\bigl{)}.$
(13)

Also, in the massless limit, the above equation reduces to Reissner–Nordström
AdS solution

$e^{2A}=e^{-2B}=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}+\frac{r^{2}}{l^{2}}.$ (14)

Now, we apply a massless limit to equation (11) and obtain charged AdS Black
Hole in $4D$ Einstein-Gauss–Bonnet gravity [19] as

$e^{2A}=e^{-2B}=1+\frac{r^{2}}{2\alpha}\Biggr{[}1\pm\sqrt{1+4\alpha\biggl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}\biggl{\\}}}\Biggr{]}.$
(15)

To find the position of the event horizon of charged AdS Black Hole in $4D$
Einstein-Gauss–Bonnet gravity, we set equation (15) equal to zero and obtain
[19]

$1-\frac{2M}{r}+\frac{Q^{2}+\alpha}{r}+\frac{r^{2}}{l^{2}}=0.$ (16)

In the absence of cosmological constant, we obtain

$r_{\pm}=M\pm\sqrt{M^{2}-Q^{2}-\alpha}.$ (17)

For Einstein-Gauss–Bonnet massive gravity with a nonvanishing cosmological
constant the expression for $r_{+}$ is complicated, so we do not represent it
here. From equation (17), we can say that the black hole solution in $4D$ EGB
massless gravity exists if and only if $M>M_{*}$ with
$M_{*}^{2}=Q^{2}+\alpha$. We will loosely follow the condition $M>M_{*}$ for
charged AdS Einstein-Gauss–Bonnet massive gravity black holes. In Fig. 1, we
plot the metric function (negative branch) of charged AdS
Einstein–Gauss–Bonnet massive gravity black holes for different values of
$\alpha$ and M. From Fig. 1 (a) and Fig. 1 (b), it is clear that the black
hole has two horizons, as the value of graviton mass increases the position of
the outer horizon increases. The position of the event horizon increases as
massive gravity parameters increase. In table 1 two roots of the metric
function (11) are estimated, the position of the horizon slowly decreases as
Gauss–Bonnet coupling parameter increases. In Fig. 1 (c) and Fig. 1 (d), we
plot the metric function, and it is clear that there are no horizon and no
black hole solutions.

Fig. 1(a)
---
$m$ |  $r_{-}$ |  $r_{+}$
0.0 |  1.1118 |  2.4489
0.5 |  1.1337 |  2.4867
1.0 |  1.2061 |  2.6419
1.5 |  1.3470 |  3.1853
Fig. 1(b)
$m$ |  $r_{-}$ |  $r_{+}$
0.0 |  1.1597 |  2.4165
0.5 |  1.1833 |  2.4532
1.0 |  1.2609 |  2.6055
1.5 |  1.4091 |  3.1500
Table 1:

(a)  $M=5$ and $\alpha=0.5$ (b) $M=5$ and $\alpha=0.8$ (c) $M=3$ and
$\alpha=0.5$ (d) $M=1$ and $\alpha=0.5$

Figure 1: $m=0.0$ denoted by solid green line, $m=0.5$ denoted by Solid black
line, $m=1.0$ denoted by dash red line and $m=1.5$ denoted by dash dot blue
line with $l=2$, $Q=3$, $c=1$, $c_{1}=-1$ and $c_{2}=1$.

In Fig. 2, we plot the metric function (negative branch) of charged AdS
Einstein–Gauss–Bonnet massive gravity in $4D$ for different values of the
charge. From the figure, it is clear that the position of the outer horizon is
the smallest for the massless case. As we increase the charge and keep the
graviton mass fixed, the position of the outer horizon decrease but the
position of the outer horizon is still greater than the massless one which is
represented by the solid green line.

Figure 2: $m=0.0$ and $Q=2$ denoted by solid green line, $m=1$ and $Q=2$
denoted by Solid black line, $m=1$ and $Q=3$ denoted by dash red line, $m=1$
and $Q=4$ denoted by dash dot blue line with $M=10$, $\alpha=0.5$ with $l=2$,
$c=1$, $c_{1}=-1$ and $c_{2}=1$.

##  III Black Hole Thermodynamics

In this section, we study the thermodynamics of charged AdS black holes in
$4D$ Einstein-Gauss–Bonnet massive gravity. The physical mass of the black
holes can be obtained from the metric function (11) by setting
$\left.e^{-2B}\right|_{r=r_{+}}=0$ as

$M=\frac{1}{2r_{+}}\biggr{[}\frac{r_{+}^{4}}{l^{2}}+r_{+}^{2}+Q^{2}+\alpha+m^{2}r_{+}^{2}\Bigl{(}\frac{cc_{1}r_{+}}{2}+c^{2}c_{2}\Bigl{)}\biggr{]}.$
(18)

The Hawking temperature of the black holes can be obtained from the relation

$T_{H}=\frac{f^{\prime}(r_{+})}{4\pi}.$ (19)

For the metric (11), Hawking temperature reads

$T_{H}=\frac{3r_{+}^{4}+l^{2}\bigl{(}r_{+}^{2}-Q^{2}-\alpha+m^{2}r_{+}^{2}(cc_{1}r_{+}^{3}+c^{2}c_{2})\bigl{)}}{4\pi
r_{+}l^{2}(r_{+}^{2}+2\alpha)}.$ (20)

In the massless limit, the Hawking temperature of charged AdS black holes in
Einstein-Gauss–Bonnet $4D$ massive gravity reduces to the Hawking temperature
[19] of charged AdS black holes in Einstein–Gauss–Bonnet $4D$ gravity as

$T_{H}=\frac{3r_{+}^{4}+l^{2}\bigl{(}r_{+}^{2}-Q^{2}-\alpha\bigl{)}}{4\pi
r_{+}l^{2}(r_{+}^{2}+2\alpha)}.$ (21)

If we take the limit $\alpha\to 0$, the Hawking temperature (21) reduces to
the Hawking temperature of Reissner–Nordström AdS black holes as

$T_{H}=\frac{3r_{+}^{4}+l^{2}\bigl{(}r_{+}^{2}-Q^{2}\bigl{)}}{4\pi
r_{+}^{3}l^{2}}.$ (22)

If we further take the chargeless limit, then the above equation reduces to
the Hawking temperature of Schwarzschild AdS black holes.

(a) $\alpha=0.5$, $Q=1$ and $c_{1}=0$ (b) $\alpha=0.8$, $Q=1$ and $c_{1}=0$

Figure 3: $m=0$ denoted by solid green line, $m=1.5$ denoted by dash black
line, $m=4.0$ denoted by dash red line, $m=5.5$ denoted by dash dot blue line
and $m=10.0$ denoted by dash dot gold line with $l=2$, $c=1$ and $c_{2}=1$.

(a) $\alpha=0.5$, $Q=5.0$ and $c_{1}=-1$ (b) $\alpha=0.5$, $Q=10.0$ and
$c_{1}=-1$

Figure 4: $m=0.0$ denoted by solid green line, $m=1.5$ denoted by dash black
line, $m=4.0$ denoted by dash red line, $m=5.5$ denoted by dash dot blue line
and $m=10.0$ denoted by dash dot gold line with $l=2$, $c=1$ and $c_{2}=1$.

In Fig. 3 and Fig. 4, we plot the Hawking temperature of the black holes with
respect to $r_{+}$ for different values of $\alpha$ and charge. In Fig. 3 (a)
and Fig. 3 (b), the Hawking temperature is plotted for different values of
$\alpha$. From the Fig., it is clear that for a critical value of horizon
radius (say, $r_{+}^{min}$) the Hawking temperature is zero, and if we
increase the horizon radius from $r_{+}^{min}$ then Hawking temperature
increases. Further increase of horizon radius from $r_{+}^{min}$ Hawking
temperature leads to attaining a local maximum for a particular value of
$r_{+}$ (say, $r_{+}^{b}$) and local maxima are slowly getting absent as
graviton mass decreases. Hawking temperature attains a minimum for a
particular value of horizon radius (say, $r_{+}^{a}$ and
$r_{+}^{a}>r_{+}^{b}>r_{+}^{min}$) and the minima slowly disappear if we
decrease the graviton mass. After attaining the minima if we further increase
the horizon radii then Hawking temperature again increases. In Fig. 4 (a) and
Fig. 4 (b), we plot the Hawking temperature with respect to the black hole
horizon for different values of charge. From Fig., it is clear that the
behavior is the same as Fig. 3 but the effects of increasing charge are local
maxima and the minima now disappear. To observe the local maxima and the
minima, we have to increase the graviton mass further, as shown by the gold
dashed dot line in Fig. 4.

In Fig. 5, the effects of charge on the Hawking temperature are shown. The
inclusion of charge slowly decreases the position of local maxima and minima.

To find the entropy of the black hole, we use the relation $dM=T_{H}dS$. Using
the Hawking temperature and mass of the black hole, we obtain the entropy as

$S=\pi r_{+}^{2}+4\pi\alpha\ln(r_{+})+S_{0},$ (23)

where $S_{0}$ is integration constant. From the above equation, one can say
that the inclusion of electric charge has no effect on the entropy of the
black hole. To derive the first law of black hole thermodynamics, we treat the
massive gravity parameters $c_{1}$ and $c_{2}$ as thermodynamics variables,
and the corresponding potential is $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$.
Apart from that, the potential corresponding to Einstein-Gauss–Bonnet
parameter $\alpha$ is $\mathcal{A}$. The thermodynamic pressure is defined as
$P=3/8\pi l^{2}$. Therefore, the first law of black hole thermodynamics in
extended phase space takes the following form:

$dM=T_{H}dS+\Phi
dQ+VdP+{\mathcal{A}}d\alpha+{\mathcal{C}_{1}}dc_{1}+{\mathcal{C}_{2}}dc_{2}.$
(24)

Now, from the first law of black hole, one can find the potential and volume
as

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\biggl{(}\frac{\partial{M}}{\partial{Q}}\biggl{)}_{S,P,\alpha,c_{1},c_{2}}=\frac{Q}{r_{+}},$
(25) $\displaystyle V$ $\displaystyle=$
$\displaystyle\biggl{(}\frac{\partial{M}}{\partial{P}}\biggl{)}_{S,Q,\alpha,c_{1},c_{2}}=\frac{4}{3}\pi
r_{+}^{3},$ (26) $\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\biggl{(}\frac{\partial{M}}{\partial{\alpha}}\biggl{)}_{S,Q,P,c_{1},c_{2}}=\frac{1}{2r_{+}},$
(27) $\displaystyle\mathcal{C}_{1}$ $\displaystyle=$
$\displaystyle\biggl{(}\frac{\partial{M}}{\partial{c_{1}}}\biggl{)}_{S,Q,P,\alpha,c_{2}}=\frac{cm^{2}r_{+}^{2}}{4},$
(28) $\displaystyle\mathcal{C}_{2}$ $\displaystyle=$
$\displaystyle\biggl{(}\frac{\partial{M}}{\partial{c_{2}}}\biggl{)}_{S,Q,P,\alpha,c_{1}}=\frac{c^{2}m^{2}r_{+}}{2}.$
(29)

##  IV Global Stability: Gibbs Free Energy

To study the global stability of the black holes we find Gibbs free energy

$G=M-T_{H}S-Q\Phi.$ (30)

Using equations (18), (20), (23) and (24), we obtain

$\displaystyle G$ $\displaystyle=$
$\displaystyle\frac{2r_{+}^{4}+l^{2}\bigl{(}2r_{+}^{2}+2Q^{2}+2\alpha+m^{2}r_{+}^{2}(cc_{1}r_{+}+2c^{2}c_{2})\bigl{)}}{4r_{+}l^{2}}$
(31) $\displaystyle-$
$\displaystyle\frac{\left(3r_{+}^{4}+l^{2}\left(r_{+}^{2}-Q^{2}-\alpha+m^{2}r_{+}^{2}(cc_{1}r_{+}+c^{2}c_{2})\right)\right)\left(4\pi\alpha\ln(r_{+})+\pi
r_{+}^{2}\right)}{4\pi r_{+}l^{2}(r_{+}^{2}+2\alpha)}-\frac{Q^{2}}{r_{+}}.$

In Fig. 5 and Fig. 6, we plot the Gibbs free energy for different values of
$\alpha$ and charge of the black hole. From Fig. 5, it is clear that the Gibbs
free energy is zero for two critical values of horizon radius (namely,
$r_{+}^{c}$, $r_{+}^{d}$ and $r_{+}^{c}>r_{+}^{d}$). The Gibbs free energy is
positive between $r_{+}^{c}$ and $r_{+}^{d}$. The positive part of the Gibbs
free energy increases as the graviton mass increases and the positive part of
the Gibbs free energy decrease as the graviton mass decrease. The positive
part of the Gibbs free energy attains its smallest value in the massless
limit. If we further increase the horizon radius $r_{+}>r_{+}^{d}$ then Gibbs
free energy goes to the negative value. In Fig. 6, we plot the effects of
charge on the Gibbs free energy. Keeping the graviton mass small $m\leq 2$, if
we increase the charge of the black hole then Gibbs free energy is completely
negative (Fig. 6(a) and 6(b)). The solid cyan points represent $r_{+}^{d}$ and
solid black points represent $r_{+}^{c}$. Finally, one can say that the
positive part of the Gibbs free energy slowly disappears due to the inclusion
of charge for small graviton mass $m\leq 2$. As graviton mass increases from
$m=2$ the behavior is similar to that in Fig. 5.

---
Figure 5: Plot of Gibbs free energy vs. horizon radius with different values
of graviton mass ($m=0,1,2,3,4$) with the fixed value of Gauss-Bonnet coupling
$(\alpha=0.1,.3,0.5,0.8)$ and $l=2$, $c=1$, $c_{1}=-1$ and $c_{2}=1$.
---
Figure 6: Plot of heat capacity vs. horizon radius with different values of
graviton mass ($m=0,1,2,3,4$) with the fixed value of Gauss-Bonnet coupling
$(\alpha=0.1,.3,0.5,0.8)$ and $l=2$, $c=1$, $c_{1}=-1$ and $c_{2}=1$.

##  V Local Stability: Heat Capacity

In this section, we study the local thermodynamical stability of the black
hole. We compute the specific heat of the black holes. The local stability of
Einstein-Gauss–Bonnet $4D$ AdS black holes is studied in Ref. [21]. The
specific heat of black holes in nonlinear electrodynamics investigated in
Refs. [32]–[34]. The local thermodynamical stability of the black holes can be
analyzed from the sign of the specific heat. If heat capacity $C_{\Phi}<0$
then the black holes are thermodynamically unstable and for $C_{\Phi}>0$ then
black holes are thermodynamically stable. The heat capacity of the black holes
is defined as

$C_{\Phi}=T_{H}\Biggl{(}\frac{dS}{dT_{H}}\Biggl{)}_{\Phi}.$ (32)

Using equations (20) and (23), this reads

$C_{\Phi}=\frac{2\pi(r_{{+}}^{2}+2\alpha)^{2}\Biggl{(}3r_{{+}}^{4}+\biggl{(}r_{{+}}^{2}-Q^{2}-\alpha+m^{2}r_{+}^{2}\bigl{(}cc_{1}r_{{+}}+c^{2}c_{2}\bigl{)}\biggl{)}l^{2}\Biggl{)}}{\Bigl{(}-(c^{2}c_{2}m^{2}+1)r_{{+}}^{4}+4cm^{2}c_{1}\alpha
r_{{+}}^{3}+(2\alpha
c^{2}c_{2}m^{2}+3Q^{2}+5\alpha)r_{{+}}^{2}+2Q^{2}\alpha+2\alpha^{2}\Bigl{)}l^{2}+3r_{{+}}^{6}+18\alpha
r_{{+}}^{4}}.$ (33)

(a) $\alpha=0.1$ (b) $\alpha=0.5$

Figure 7: $m=0.0$ denoted by solid green line, $m=1.0$ denoted by dash black
line, $m=2.0$ denoted by dash red line, $m=3.0$ denoted by dash dot blue line
and $m=4.0$ denoted by dash dot gold line with $Q=1$, $l=3$, $c=1$, $c_{1}=-1$
and $c_{2}=1$.

In Fig. 7, we plot the specific heat for two different values of the Einstein-
Gauss-Bonnet coupling parameter and it is discontinuous at some critical value
of horizon radius for a larger value of graviton mass ($m>1$) which indicates
that a second-order phase transition occurs for charged AdS black holes in
$4D$ Einstein-Gauss–Bonnet massive gravity. For the smaller value of graviton
mass($m=1$), no divergences occur but the smaller-size black holes are
thermodynamically unstable as specific heat is negative. As the size of the
black hole increases a phase transition occurs, i.e., the specific heat of the
black hole changes from a negative value to a positive value. If we take
graviton mass as zero then similar kinds of phenomena occur.

In Ref. [50], the specific heat of EGB massive gravity black hole is studied
with $Q=0$. In the chargeless case, $Q=0$ two diverging points appear for two
critical values of horizon radius which separate three regions, i.e. two-
second order phase transition occurs for such a black hole. Between two
diverging points, specific heat is negative, which indicates that the black
hole is thermodynamically unstable in this region. The inclusion of charge
removes one diverging point and we are left with only one diverging point,
i.e. in the case of charged black hole only one second-order phase transition
occurs. For $m=0$ behavior of specific heat same as $Q=0$ and $Q\neq 0$.

##  VI Dynamic Stability: Quasinormal modes

One of the methods to study the dynamic stability of the black holes is
studying the nature of quasinormal modes (QNM) which are characterized by
complex numbers. If the imaginary part of QNM is positive, the black hole is
unstable; however, if negative, the black hole is stable.

We compute the QNM and quasinormal frequency (QNF) of the above black hole
using the scalar field perturbation. We need to consider the scalar field
$\Phi$ in the background of the black hole (11). The equation of these
perturbations takes the form

$\frac{1}{\sqrt{-g}}\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\right)\Phi=0,$
(34)

where $g^{\mu\nu}$ are the metric components of the metric (4). The mode
decomposition of the scalar perturbation in terms of spherical harmonics is
given by

$\Phi=\frac{1}{r}\sum_{lm}e^{i\omega t}\phi_{lm}Y^{m}_{l}(\theta,\phi),$ (35)

here $l$, $m$, $Y^{m}_{l}$, and $\omega$ are respectively the angular quantum
number, magnetic quantum numbers, spherical harmonic, and the oscillating
frequency of the scalar field. Substituting the value of $\Phi$ in Eq. (34)
and using the tortoise coordinate $dr^{*}=dr/e^{2A}$ , we get the Schrodinger-
like form,

$\left(\frac{d^{2}}{dr^{*^{2}}}+\omega^{2}-V(r^{*})\right)\phi=0,$ (36)

where $V(r^{*})$ is the effective potential and has the form

$V(r^{*})=e^{2A}\left({A^{\prime}(r)}+\frac{l(l+1)}{r^{2}}\right),$ (37)

where $l$ is the harmonic index. To find the QNF one has to impose boundary
conditions near the event horizon. These boundary conditions can be written as

$\displaystyle\phi(r_{\star})\to e^{i\omega r_{\star}},\qquad\qquad
r_{\star}\to-\infty,$ (38) $\displaystyle\phi(r_{\star})\to e^{-i\omega
r_{\star}},\qquad\qquad r_{\star}\to\infty,$ (39)

where the $+$ sign corresponds to ingoing waves at the horizon and $-$ sign
corresponds to outgoing waves at the infinity. The frequencies corresponding
to the QNM are given by $\omega=\omega_{R}+i\omega_{I}$, whose $\omega_{R}$
and $\omega_{I}$ are the oscillating damping components of the frequency. We
use the WKB approximation to find the QNMs and QNFs of the obtained black hole
solution (11). The WKB formula has the form [67, 68, 69]

$i\frac{\omega^{2}-V_{0}}{\sqrt{-2V^{\prime\prime}_{0}}}=n+\frac{1}{2}.$ (40)

where $V_{0}$ is the height of the barrier and $V^{\prime\prime}_{0}$ is the
second derivative of the potential with respect to the tortoise coordinate.
The numerical value of QNM and QNF for different values of graviton mass is
depicted in Tab. 2

| $\alpha=0.1$ | $\alpha=0.2$ | $\alpha=0.3$
---|---|---|---
$m$ | $\omega=\omega_{R}+i\omega_{I}$ | $\omega=\omega_{R}+i\omega_{I}$ | $\omega=\omega_{R}+i\omega_{I}$
1 | 0.0036 - 0.0036 $i$ | 0.0011 - 0.001 $i$ | 0.0012 - 0.0012 $i$
2 | 0.0082 - 0.0091 $i$ | 0.0036 - 0.0039 $i$ | 0.0036 - 0.0038 $i$
3 | 0.0065 - 0.0074 $i$ | 0.0042 - 0.0047 $i$ | 0.0054 - 0.0058 $i$
4 | 0.0052 - 0.0060 $i$ | 0.0056 - 0.0062 $i$ | 0.0034 - 0.0036 $i$
5 | 0.0065 - 0.0071 $i$ | 0.0043 - 0.0046 $i$ | 0.0052 - 0.0055 $i$
Table 2: The numerical values of QNMs with different values of gravitation
mass $(m)$ and Gauss-Bonnet coupling $(\alpha)$ with a fixed value of
$M=1,Q=1,c=1,c_{1}=-1,c_{2}=1,n=1$, and $l=10$.

From the table 2, we can see clearly that the imaginary part of the QNMs in
the obtained black hole solution (11) is negative. So the black hole solution
is stable.

##  VII Van der Waals Like Phase Transition

In this section, we study the phase transition of charged AdS black holes in
$4D$ Einstein-Gauss–Bonnet massive gravity. The phase transition of the black
holes in massive gravity is studied in Refs. [70, 50]. Hawking-Page phase
transition for static and rotating $4D$ Gauss–Bonnet black hole is studied in
Refs. [71, 22]. The phase transition of charged AdS black hole in
Einstein–Gauss–Bonnet massless gravity is studied in Ref. [21]. Phase
transition of charged AdS $4D$ Einstein-Gauss–Bonnet black hole in nonlinear
electrodynamics is also studied [32]. The Van der Waals equation of state for
real fluids is given by

$P=\frac{T}{v-b}-\frac{a}{v^{2}},$ (41)

where $v$ is the specific volume of the fluids, $a$ represents the interaction
between the molecules of the fluids and $b$ describes the non-zero of the
molecules. Now, from the Hawking temperature (20), we obtain

$P=\frac{T}{v}+\frac{8T\alpha}{v^{3}}+\frac{2Q^{2}}{v^{4}\pi}-\frac{1}{2v^{2}\pi}+\frac{2\alpha}{v^{4}\pi}-\frac{c^{2}c_{2}m^{2}}{2v^{2}\pi}-\frac{cc_{1}m^{2}}{4v\pi},$
(42)

where specific volume $v$ is defined by [72]

$v=\frac{6V}{A}\approx 2r_{+},$ (43)

where $A$ is the area of the black hole. To obtain the critical points, we use
the following conditions:

$\Biggl{(}\frac{\partial{P}}{\partial{v}}\Biggl{)}_{T_{c},v_{c}}=\Biggl{(}\frac{\partial^{2}{P}}{\partial{v}^{2}}\Biggl{)}_{T_{c},v_{c}}=0.$
(44)

Using equations (42) and (44), we obtain the condition for critical volume as

$-2(c^{2}c_{2}m^{2}+1)v_{c}^{4}+24\alpha cc_{1}m^{2}v_{c}^{3}+48(\alpha
c^{2}c_{2}m^{2}+Q^{2}+2\alpha)v_{c}^{2}+384Q^{2}\alpha+384\alpha^{2}=0.$ (45)

Equation (45) can not be solved analytically, we numerically solved the above
equation and estimate the critical points as shown in the tables below.

m |  ${v_{c}}$ |  $P_{c}$ |  $T_{c}$ |  $\rho_{c}$
---|---|---|---|---
0.0 |  5.4328 |  0.0025 |  0.0380 |  0.3702
0.1 |  5.4030 |  0.0026 |  0.0369 |  0.3806
0.2 |  5.3161 |  0.0027 |  0.0337 |  0.4259
0.3 |  5.1788 |  0.0030 |  0.0285 |  0.5451
0.4 |  5.0010 |  0.0034 |  0.0213 |  0.7982
0.5 |  4.7940 |  0.0040 |  0.0123 |  1.5590
Table 3: Values of critical volume ($v_{c}$), critical pressure ($P_{c}$), critical temperature ($T_{c}$) and $\rho_{c}=P_{c}v_{c}/T_{c}$ for different graviton mass with $Q=1$, $\alpha=0.1$, $c=1$ $c_{1}=-2$ and $c_{2}=0.75$. Q |  ${v_{c}}$ |  $P_{c}$ |  $T_{c}$ |  $\rho_{c}$
---|---|---|---|---
0.0 |  1.9620 |  0.0206 |  0.0733 |  0.5513
0.2 |  2.1519 |  0.0177 |  0.0659 |  0.5779
0.4 |  2.6326 |  0.0125 |  0.0501 |  0.6568
0.6 |  3.2763 |  0.0083 |  0.0345 |  0.7882
0.8 |  4.0092 |  0.0057 |  0.0219 |  1.0434
1.0 |  4.7940 |  0.0040 |  0.0123 |  1.5577
Table 4: Values of critical volume($v_{c}$), critical pressure ($P_{c}$), critical temperature ($T_{c}$) and $\rho_{c}=P_{c}v_{c}/T_{c}$ for different charged of the black hole with $m=0.5$, $\alpha=0.1$, $c=1$ $c_{1}=-2$ and $c_{2}=0.75$. $\alpha$ |  ${v_{c}}$ |  $P_{c}$ |  $T_{c}$ |  $\rho_{c}$
---|---|---|---|---
0.0 |  0.4827 |  0.3516 |  0.4464 |  0.3801
0.1 |  2.2693 |  0.0132 |  0.0778 |  0.3834
0.2 |  3.1460 |  0.0068 |  0.0543 |  0.3939
0.3 |  3.8142 |  0.0046 |  0.0437 |  0.4014
0.4 |  4.3728 |  0.0035 |  0.0373 |  0.4103
0.5 |  4.8610 |  0.0028 |  0.0330 |  0.4124
Table 5: Values of critical volume ($v_{c}$), critical pressure ($P_{c}$),
critical temperature ($T_{c}$) and $\rho_{c}=P_{c}v_{c}/T_{c}$ for different
Gauss–Bonnet coupling parameter with $Q=0.1$, $m=0.2$, $c=1$ $c_{1}=-2$ and
$c_{2}=0.75$.

In tables 3, 4 and 5, we numerically solve equation (45) for different values
of graviton mass, charge and Einstein-Gauss-Bonnet coupling parameter, we
estimate the value of critical volume ($v_{c}$), critical pressure ($P_{c}$),
critical temperature ($T_{c}$) and $\rho_{c}$. From table 3, we can say that
as graviton mass increases from zero then the critical volume ($v_{c}$),
critical temperature ($T_{c}$) decreases, critical pressure ($P_{c}$), and
$\rho_{c}$ increases. The effects of black hole charge on the critical
parameters are shown in table 4, keeping the graviton mass fixed. As the
charge of the black holes increases critical volume ($v_{c}$) and $\rho_{c}$
increase, however, critical pressure and critical temperature decrease. The
effects of the Gauss–Bonnet coupling parameter ($\alpha$) on the critical
parameters are shown in table 5. In Fig. 7, we plot the Hawking temperature
for different values of Gauss–Bonnet coupling parameter and charge with
$P<P_{c}$, $P=P_{c}$ and $P>P_{c}$.

In Fig. 8(a) and 8(b), the Hawking temperature is depicted for different
values of $\alpha$ keeping pressure fixed. When pressure is less than critical
pressure ($P_{c}$), the curve has two critical points (one of them is maxima
and another is minima). For the pressure equal to the critical pressure
($P_{c}$), two turning points come to an inflection point and when $P>P_{c}$
the curve does not attain any turning points. The effects of charge on the
Hawking temperature are shown in Fig. 8(c) and 8(d). The inclusions of charge
basically reduced the position of local maxima and minima when $P<P_{c}$.
Furthermore, if we increase the charge then the position of local maxima and
minima decrease by a significant amount (Fig. 8(d)). The rest of the behaviour
is similar to Fig. 8(a) and Fig. 8(b). The behaviour of Hawking temperature of
charged less black hole for $P\leq P_{c}$ and $P>P_{c}$ is shown in Ref. [50].
The Hawking temperature of the black hole with $Q=0$ attains local maxima and
minima when $P<P_{c}$. In the chargeless case position of local maxima and
minima is higher than the charged black hole.

(a) $\alpha=0.5$, $Q=1$ and $c_{1}=0$ (b) $\alpha=0.8$, $Q=1$ and $c_{1}=0$
(c) $\alpha=0.5$, $Q=3.0$ and $c_{1}=-0.1$ (d) $\alpha=0.5$, $Q=5.0$and
$c_{1}=-0.1$

Figure 8: $P=0.25P_{c}$ denoted by solid black line, $P=P_{c}$ denoted by dash
red line and $P=3P_{c}$ denoted by dash dot blue line with $c=1$, $c_{2}=1$
and $m=1.5$.

In Fig. 9, we plot Gibbs free energy Vs. temperature for different values of
Gauss–Bonnet coupling parameter and black holes charge with $P<P_{c}$,
$P=P_{c}$ and $P>P_{c}$. In Fig. 9(a) and 9(b), Gibbs free energy for
different values of Gauss–Bonnet coupling parameters is depicted. When
pressure is less than the critical pressure ($P_{c}$), Gibbs free energy shows
swallow tail (triangular shape) behavior, which indicates that the system
undergoes a first-order phase transition, i.e., below the critical pressure a
transition between small black hole (SBH) and large black hole (LBH) occurs.
The Gibbs free energy of LBH is smaller compared to the SBH. At the point of
intersection of the curve ($P<P_{c}$), where first-order phase transition
occurs, the entropy of the system is discontinuous as entropy depends on the
horizon radius of the black holes and the radius of the SBH and LBH is
different. For $P=P_{c}$, the swallow tail behavior disappears at which a
second-order phase transition occurs. For $P>P_{c}$ swallow tail behavior
completely disappears and no phase transition occurs. A similar kind of
behavior is shown in Fig.9(c) and Fig.9(d) for different values of charge.

(a) $\alpha=0.5$ and $Q=1.0$ (b) $\alpha=0.8$ and $Q=1.0$ (c) $\alpha=0.5$ and
$Q=3.0$ (d) $\alpha=0.5$ and $Q=5.0$

Figure 9: $P=0.25P_{c}$ denoted by solid black line, $P=P_{c}$ denoted by dash
red line and $P=3P_{c}$ denoted by dash dot blue line with $c=1$,
$c_{1}=-0.1$, $c_{2}=1$ and $m=1$.

##  VIII Joule-Thomson Expansion

In this section, we discuss the effects of massive gravity on the Joule-
Thomson expansion of charged AdS black holes in $4D$ Einstein–Gauss–Bonnet
massive gravity. The Joule-Thomson expansion of charged AdS black holes was
first studied in Ref. [59]. After that Joule-Thomson expansion of $D$
dimensional black holes was studied in Ref. [73]. Using numerical
investigation, Joule-Thomson expansion of Kerr-AdS and Kerr-Newman-AdS is also
studied [60, 74]. Joule-Thomson expansion of charged AdS black hole in $4D$
Einstein massive gravity discussed in Ref. [66]. Joule-Thomson expansion of
charged AdS $4D$ Einstein massive gravity black hole in Maxwell and Born-
Infeld electrodynamics discussed in Ref. [21, 75]. The Joule-Thomson
thermodynamic coefficient is given by

$\mu_{J}=\biggl{(}\frac{\partial{T}}{\partial{P}}\biggl{)}_{M}=\frac{1}{C_{P}}\biggr{[}T\Bigl{(}\frac{\partial{V}}{\partial{P}}\Bigl{)}_{P}-V\biggr{]}=\frac{(\partial{T}/\partial{r_{+}})_{M}}{(\partial{P}/\partial{r_{+}})_{M}}.$
(46)

The Joule-Thomson effect is an isenthalpic process, which means that enthalpy
remains constant during the process. In the Joule-Thomson process pressure
always decrease but the temperature can increase/decrease, thus Joule-Thomson
thermodynamic coefficient ($\mu$) can be negative/positive. When $\mu>0$ the
Joule–Thomson expansion corresponds to the cooling region of the isenthalpic
or constant mass curve and $\mu<0$ corresponds to the heating region of the
isenthalpic or constant mass curve. The Joule - Thomson thermodynamic
coefficient vanishes for some particular value of temperature, which is known
as inverse temperature ($T_{i}$) and corresponding pressure is known as
inverse pressure ($P_{i}$). The cooling and heating regions are separated by
the set of points ($P_{i}$, $T_{i}$) and the curved formed by the set of
points ($P_{i}$, $T_{i}$) known as the inverse curve. Clearly, the region
above the inverse curve is known as the cooling region and the region below
the inverse curve is known as the heating region. At inverse temperature sign
of Joule–Thomson coefficient changes $\mu_{J}(T_{i})=0$. From the above
equation, we obtain inverse temperature

$T_{i}=V\biggl{(}\frac{\partial{T}}{\partial{V}}\biggl{)}_{P}=\frac{r_{+}}{3}\biggl{(}\frac{\partial{T}}{\partial{r_{+}}}\biggl{)}_{P}.$
(47)

From equation (18) and using $P=3/8\pi l^{2}$, we obtain pressure in terms of
black hole mass

$P=\frac{3}{4\pi
r_{+}^{2}}\left[\frac{M}{r_{+}}-\frac{Q^{2}}{2r_{+}^{2}}-\frac{1}{2}-\frac{\alpha}{2r_{+}^{2}}-\frac{m^{2}c^{2}c_{2}}{2}-\frac{m^{2}cc_{1}r_{+}}{4}\right].$
(48)

Now, from Hawking temperature and using the relations $P=1/8\pi l^{2}$, we
obtain

$T=\frac{1}{4\pi r(r_{+}^{2}+2\alpha)}\left[8\pi
Pr_{+}^{4}-Q^{2}+r_{+}^{2}-\alpha+c^{2}c_{2}m^{2}r_{+}^{2}+m^{2}cc_{1}r_{+}^{3}\right].$
(49)

(a) $\alpha=0.5$ (b) $\alpha=0.8$

Figure 10:  $M=2$ denoted by solid black line, $M=2.5$ denoted by red dash
line, $M=3$ denoted by blue dash-dot line with $m=1$, $Q=1$, $c=1$, $c_{1}=-1$
and $c_{2}=1$. A solid gold line denotes an inverse curve.

(a) $\alpha=0.5$ (b) $\alpha=0.8$

Figure 11:  $m=0$ denoted by solid black line, $m=1$ denoted by red dash line,
$m=2$ denoted by blue dash dot line and $m=3$ denoted by green dash dot line
with $M=5$, $Q=1$, $c=1$, $c_{1}=-1$ and $c_{2}=1$.

(a) $Q=2.0$ (b) $Q=4.0$

Figure 12: $M=20$, $\alpha=0.5$ and $c=1$. Left panel: $m=0$ denoted by the
solid black line, $m=3$ denoted by a red dash line, $m=5$ denoted by a blue
dash-dot line and $m=7$ denoted by an orange dash-dot line. Right panel: $m=0$
denoted by the solid black line, $m=1$ denoted by a red dash line, $m=2$
denoted by a blue dash-dot line and $m=3$ denoted by an orange dash-dot line.

(a) $c_{1}=0$ (b) $c_{2}=0$

Figure 13:  $M=2$, $Q=1$, $\alpha=0.5$ and $c=1$. Left panel: $c_{2}=0.1$
denoted by an orange dash line, $c_{2}=1$ denoted by a red dash line,
$c_{2}=-1$ denoted by a green dash-dot line, and $c_{2}=-2$ denoted by a blue
dash-dot line. Right panel: $c_{1}=1$ denoted by a red dash line, $c_{1}=2$
denoted by an orange dash line, $c_{1}=-1$ denoted by a green dash-dot line
and $c_{2}=-2$ denoted by a blue dash-dot line.

Using equation (48) and equation (21) constant mass curve can be obtained. In
Fig. 10-Fig. 14 we plot the constant mass curve. In Fig. 10, constant mass and
an inverse curve are shown for different values of black hole mass. The left
region of the inverse curve represents cooling and the right region represents
heating. In Fig. 11 and Fig. 12 constant mass curve is shown for different
values of Gauss–Bonnet coupling parameter and charge of the black hole. The
effects of parameters $c_{1}$ and $c_{2}$ are shown in Fig. 13 \- Fig. 15.

(a) $c_{1}=1$ and $c_{2}=1$ (b) $c_{1}=1$ and $c_{2}=-1$

Figure 14: $m=0$ denoted by solid black line, $m=0.5$ denoted by dash red
line, $m=1$ denoted by dash orange line, $m=1.5$ denoted by blue dash dot line
and $m=2$ denoted by dash dot green line with $M=2$, $Q=1$, $\alpha=0.5$ and
$c=1$.

(a) $c_{1}=-1$ and $c_{2}=1$ (b) $c_{1}=-1$ and $c_{2}=-1$

Figure 15: $m=0$ denoted by solid black line, $m=0.5$ denoted by dash red
line, $m=1$ denoted by dash orange line, $m=1.5$ denoted by blue dash dot line
and $m=2$ denoted by dash dot green line with $M=2$, $Q=1$, $\alpha=0.5$ and
$c=1$.

(a) $\alpha=0.5$ (b) $m=1$

Figure 16:  $Q=1$, $c=1$, $c_{1}=-1$ and $c_{2}=1$. Left panel: $m=0$ denoted
by the solid black line, $m=1$ denoted by red dash line, $m=2$ denoted by
orange dash line, $m=3$ denoted by blue dash-dot line, $m=4$ denoted by green
dash-dot line and $m=5$ denoted by a gold dash-dot line. Right panel:
$\alpha=0$ denoted by a solid black line, $\alpha=0.1$ denoted by red dash
line, $\alpha=0.2$ denoted by orange dash line, $\alpha=0.3$ denoted by blue
dash-dot line, $\alpha=0.4$ denoted by green dash-dot line and $\alpha=0.5$
denoted by a gold dash-dot line.

We use equations (47) and (49) and obtain the inverse pressure as

$P_{i}=\frac{6Q^{2}r_{+}^{2}-4r^{4}+8Q^{2}\alpha+2\alpha
r_{+}^{2}+8\alpha^{2}-4c^{2}c_{2}m^{2}r_{+}^{4}-3cc_{1}m^{2}r_{+}^{5}-4\alpha
c^{2}c_{2}m^{2}r_{+}^{2}-2\alpha cc_{1}m^{2}r_{+}^{3}}{16\pi r_{+}^{6}}.$ (50)

Now, using the relations (49) and (50), we obtain

$T_{i}=\frac{-m^{2}cc_{1}r_{+}^{3}+(-2c^{2}c_{2}m^{2}-2)r_{+}^{2}+4Q^{2}+4\alpha}{8\pi
r_{+}^{3}}.$ (51)

(a) $c_{2}=0$ (b) $c_{1}=0$

Figure 17: $Q=1$, $c=1$, $\alpha=0.5$ and $m=1$. Left panel: $c_{1}=2$ denoted
by a red dash line, $c_{1}=1$ denoted by an orange dash line, $c_{1}=-1$
denoted by a blue dash-dot line, $c_{1}=-2$ denoted by a green dash-dot line.
Right panel: $c_{2}=2$ denoted by a red dash line, $c_{2}=1$ denoted by an
orange dash line, $c_{2}=-1$ denoted by a blue dash-dot line, $c_{2}=-2$
denoted by a green dash-dot line.

(a) $c_{1}=0$ (b) $c_{2}=0$

Figure 18: $M=2.5$, $Q=1$, $c=1$, $\alpha=0.5$ and $m=1$. Left panel:
$c_{2}=-2$ denoted by a red dash line, $c_{2}=-1$ denoted by an orange dash
line, $c_{2}=1$ denoted by a blue dash-dot line, $c_{2}=2$ denoted by a green
dash-dot line. Right panel: $c_{1}=-2$ denoted by a red dash line, $c_{1}=-1$
denoted by an orange dash line, $c_{1}=1$ denoted by a blue dash-dot line,
$c_{1}=2$ denoted by a green dash-dot line.

Finally, we will drive the Joule - Thomson thermodynamic coefficient using
equation (46), equation (48) and equation (49)

$\mu_{J}=\frac{8r_{+}^{3}\Bigl{(}6Mr_{+}^{3}-6Q^{2}r_{+}^{2}-r_{+}^{4}-4Q^{2}\alpha-4\alpha
r_{+}^{2}-4\alpha^{2}-c^{2}c_{2}m^{2}r_{+}^{4}+2\alpha
c^{2}c_{2}m^{2}r_{+}^{2}+\alpha
cc_{1}m^{2}r_{+}^{3}\Bigl{)}}{3(r_{+}^{2}+2\alpha)^{2}\Bigl{(}12Mr_{+}-8Q^{2}-4r_{+}^{2}-8\alpha-4c^{2}c_{2}m^{2}r_{+}^{2}-m^{2}cc_{1}r_{+}^{3}\Bigl{)}}.$
(52)

(a) $c_{1}=1$ and $c_{2}=1$ (b) $c_{1}=1$ and $c_{2}=-1$ (c) $c_{1}=-1$ and
$c_{2}=1$ (d) $c_{1}=-1$ and $c_{2}=-1$

Figure 19: $m=0$ denoted by solid black line, $m=0.5$ denoted by red dash
line, $m=1$ denoted by orange dash line, $m=1.5$ denoted by blue dash-dot
line, $m=2$ denoted by green dash-dot line with $M=2.5$, $Q=1$, $c=1$, and
$\alpha=0.5$.

In Fig. 19 and Fig. 18, Joule - Thomson thermodynamic coefficient is plotted.
$\mu<0$ represent heating phase and $\mu>0$ represent cooling phase.

##  IX Conclusions

In this work, we have found an exact solution of Einstein-Gauss–Bonnet massive
gravity with charge in $4D$ AdS space and the horizon structure of the black
holes is discussed. The physical mass (enthalpy) and Hawking temperature of
the black holes were computed. Treating the cosmological constant as the
pressure, we drive the first law of black hole thermodynamics. To check global
stability, the Gibbs free energy of the black hole was computed. For the local
stability of the black holes we estimated the specific heat. Here, we found
that as the mass of the graviton increases one divergent point appears. For
$m=0,1$, the specific heat of the black hole goes from a negative value
(unstable phase) to a positive one (stable phase) and when $m>1$ a second-
order phase transition occurs. Furthermore, we investigated the Van der Walls-
like phase transition of the black holes. The effects of graviton mass, the
charge of the black holes and the Gauss–Bonnet coupling parameter on the
critical points were also studied. As the mass of the graviton increase,
critical pressure also increases and critical temperature shows the opposite
behavior. When charge and Gauss–Bonnet coupling parameters of the black hole
increase, we found that both critical pressure and temperature decrease. From
the $G-T$ plot, we observed that a swallow tail appears below the critical
point which indicates a first-order phase transition, and at the critical
point, a second-order phase transition occurs. Finally, Joule-Thomson
expansion of charged AdS black hole in $4D$ Einstein–Gauss–Bonnet massive
gravity was studied. The effects of the Gauss–Bonnet coupling parameter and
massive gravity parameters on the constant mass and inverse curve are also
depicted.

## Acknowledgement

D.V.S. thanks University Grant Commission for the Start-Up Grant No.
30-600/2021(BSR)/1630.

## Data Availability Statement

No Data is associated with the manuscript.

## References

  * [1] C. Barrabes and G. F. Bressange, Class. Quant. Grav. 14, 805 (1997).
  * [2] R.-G. Cai and Y. S. Myung, Phys. Rev. D 56, 3466 (1997).
  * [3] S. Capozziello and A. Troisi, Phys. Rev. D 72, 044022 (2005).
  * [4] T. P. Sotiriou, Class. Quant. Grav. 23, 5117 (2006).
  * [5] J. W. Moffat, JCAP 03, 004 (2006).
  * [6] V. Faraoni, Phys. Rev. D 75, 067302 (2007).
  * [7] D. Lovelock, J. Math. Phys. 12, 498 (1971).
  * [8] D. Lovelock, J. Math. Phys. 13, 874 (1972).
  * [9] N. Deruelle and L. Farina-Busto, Phys. Rev. D 41, 3696 (1990).
  * [10] T. Tangphati, A. Pradhan, A. Banerjee and G. Panotopoulos, Physics of the Dark Universe 33, 100877 2021).
  * [11] J.M.Z. Pretel, A. Banerjee and A. Pradhan, Eur. Phys. J. C 82, 180 (2022).
  * [12] T. Tangphati, A. Pradhan, E. Abdelghani and A. Banerjee, Physics Letters B 819, 136423 (2021).
  * [13] T. Tangphati, A. Pradhan, E. Abdelghani and A. Banerjee, Physics Letters B 819, 136423 (2021).
  * [14] J. M. Cline and H. Firouzjahi, Phys. Rev. D 64, 023505 (2001).
  * [15] T. Nihei, N. Okada and O. Seto, Phys. Rev. D 71, 063535 (2005).
  * [16] M. Demetrian, Gen. Rel. Grav. 38, 953 (2006).
  * [17] C. Lanczos, Annals Math. 39, 842 (1938).
  * [18] D. Glavan and C. Lin, Phys. Rev. Lett. 124, 081301 (2020).
  * [19] P. G. S. Fernandes, Phys. Lett. B 805, 135468 (2020).
  * [20] P. G. S. Fernandes and P. Carrilho, T. Clifton and D. J. Mulryne, Class. Quant. Grav. 39, 063001 (2022).
  * [21] K. Hegde, A. Naveena Kumara, C. L. Ahmed Rizwan, Ajith K.M. and Md S.Ali, 2003.08778.
  * [22] Y. Y. Wang, B.-Yu Su and N. Li, Phys. Dark Univ. 31, 100769 (2021).
  * [23] S.-W. Wei and Y.-X. Liu, Phys. Rev. D 101, 104018 (2020).
  * [24] D. V. Singh, S. G. Ghosh and S. D. Maharaj, Phys. Dark Univ. 30, 100730 (2020).
  * [25] D. V. Singh and S. Siwach, Phys. Lett. B 808, 135658 (2020).
  * [26] B. Eslam Panah, K. Jafarzade and S. H. Hendi, Nucl. Phys. B 961, 115269 (2020).
  * [27] D. V. Singh, B. K. Singh and S. Upadhyay, Annals Phys. 434, 168642 (2021).
  * [28] N. Godani, D. V. Singh and G. C. Samanta, Phys. Dark Univ. 35, 100952 (2022).
  * [29] S. H. Hendi and S. Panahiyan, Phys. Rev. D 90, 124008 (2014).
  * [30] K. Yang and B.-M. Gu, S.-W. Wei and Y.-X. Liu, Eur. Phys. J. C 80, 662 (2020).
  * [31] D. V. Singh, V. K. Bhardwaj and S. Upadhyay, Eur. Phys. J. Plus 137, 969 (2022).
  * [32] S. G. Ghosh, D. V. Singh, R. Kumar and S. D. Maharaj, Annals Phys. 424, 168347 (2021).
  * [33] S. I. Kruglov, Symmetry 13, 204 (2021).
  * [34] S. I. Kruglov, Annals Phys. 428, 168449 (2021).
  * [35] S. I. Kruglov, EPL 133, 6 (2021).
  * [36] S. I. Kruglov, Symmetry 13, 944 (2021).
  * [37] S. I. Kruglov, Universe 7, 249 (2021).
  * [38] B. P. Abbott et al., Phys. Rev. Lett. 116, 221101 (2016).
  * [39] M. Fierz and W. Pauli, Proc. Roy. Soc. Lond. A 173, 211 (1939).
  * [40] M. Fierz, Helv. Phys. Acta 12, 3 (1939).
  * [41] D. G. Boulware and S. Deser, Phys. Rev. D 6, 3368 (1972).
  * [42] C. de Rham and G. Gabadadze, Phys. Rev. D 82, 044020 (2010).
  * [43] C. de Rham, G. Gabadadze and A. J. Tolley, Phys. Rev. Lett. 106, 231101 (2011).
  * [44] S. H. Hendi, S. Panahiyan and B. Eslam Panah, JHEP 01, 129 (2016).
  * [45] S. H. Hendi, B. Eslam Panah and S. Panahiyan, Class. Quant. Grav. 33, 235007 (2016).
  * [46] S. H. Hendi and M. Momennia, JHEP 10, 207 (2019).
  * [47] B. K. Singh and R. P. Singh and D. V. Singh, Eur. Phys. J. Plus 135, 862 (2020).
  * [48] S. H. Hendi, S. Panahiyan, S. Upadhyay and B. Eslam Panah, Phys. Rev. D 95, 084036 (2017).
  * [49] S. Upadhyay, S. H. Hendi, S. Hossein, S. Panahiyan and B. Eslam Panah, PTEP 2018, 093E01 (2018).
  * [50] S. Upadhyay, and D. V. Singh, Eur. Phys. J. Plus 137, 383 (2022).
  * [51] S W. Hawking and D. N. Page, Commun. Math. Phys. 87, 577 (1983).
  * [52] S. Wang, S.-Q. Wu, F. Xie and L. Dan, Chin. Phys. Lett. 23, 1096- (2006).
  * [53] D. Kastor, S. Ray and J. Traschen, Class. Quant. Grav. 26, 195011 (2009).
  * [54] D. Kastor, S. Ray and J. Traschen, Class. Quant. Grav. 27, 235014 (2010).
  * [55] B. P. Dolan, Class. Quant. Grav. 28, 125020 (2011).
  * [56] B. P. Dolan, Class. Quant. Grav. 28, 235017 (2011).
  * [57] D. Kubizňák and R. B. Mann, JHEP 07, 033(2012).
  * [58] S. Gunasekaran and R. B. Mann D. Kubiznak, JHEP 11, 110 (2012).
  * [59] O. Ökcü and E. Aydıner, Eur. Phys. J. C 77, 24(2017).
  * [60] O. Ökcü and E. Aydıner, Eur. Phys. J. C 78, 123 (2018).
  * [61] G. G. L. Nashed and S. Capozziello, Eur. Phys. J. C 79, 911 (2019).
  * [62] G. G. L. Nashed and S. Capozziello, Int. J. Mod. Phys. A 33, 1850076 (2018).
  * [63] G. G. L. Nashed and S. Capozziello, Phys. Rev. D 99, 104018 (2019).
  * [64] C. H. Nam, Eur. Phys. J. Plus 135, 259 (2020).
  * [65] R.-G. Cai, Y.-P. Hu, Q.-Y. Pan and Y.-L. Zhang, Phys. Rev. D 91, 024032 (2015).
  * [66] C. H. Nam, Eur. Phys. J. Plus 135, 259 (2020).
  * [67] B. F. Schutz and C. M. Will, The Astrophysical Journal 291, L33-L36, (1985).
  * [68] S. Iyer and C. M. Will, Phys. Rev. D 35, 3621 (1987).
  * [69] R. A. Konoplya, Phys. Rev. D 68, 024018 (2003).
  * [70] S. Fernando, Mod. Phys. Lett. A 31, 1650096 (2016).
  * [71] B.-Yu Su, Y. Y. Wang and N. Li, Eur. Phys. J. C 80, 305 (2020).
  * [72] A. Rajagopal, D. Kubizňák and R.B. Mann, Phys. Lett. B 737, 277 (2014).
  * [73] J.-X. Mo, G.-Q. Li, S.-Q. Lan and X.-B. Xu, Phys. Rev. D 98, 124032 (2018).
  * [74] Z.-W. Zhao, Y.-H. Xiu, and N. Li, Phys. Rev. D 98, 124003 (2018).
  * [75] C.-M. Zhang, M. Zhang and D.-C. Zou, Mod. Phys. Lett. A 37, 2250063 (2022).

## Appendix A

The Ricci and Kretschmann scalar for the metric function (11) are given by

$\displaystyle R$ $\displaystyle=$
$\displaystyle-\frac{f_{{rr}}(r)r^{2}+4f_{r}(r)r+2f(r)-2}{r^{2}},$ (53)
$\displaystyle K$ $\displaystyle=$ $\displaystyle
f_{{rr}}^{2}(r)+\frac{4f_{r}^{2}(r)}{r^{2}}+\frac{4(-1+f(r))^{2}}{r^{4}},$
(54)

where we use $f(r)=e^{2A(r)}$, $f_{rr}(r)$ stands for $d^{2}f(r)/dr^{2}$ and
$f_{r}(r)$ stands for $df(r)/dr$.

The differentiation of $f(r)$ with respect to $r$ is given by

$\displaystyle f_{r}$ $\displaystyle=$
$\displaystyle\frac{r}{\alpha}\left[1-\sqrt{1+4\alpha\left\\{\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\right\\}}\right]$
(55) $\displaystyle-$
$\displaystyle\frac{r^{2}\Bigr{[}-\frac{6M}{r^{4}}+\frac{4Q^{2}}{r^{5}}-\frac{m^{2}cc_{1}}{2r^{2}}+\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{r^{3}}\Bigr{]}}{\sqrt{1+4\alpha\biggl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\biggl{\\}}}},$
$\displaystyle f_{rr}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha}\Biggr{[}1-\sqrt{1+4\alpha\biggl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\biggl{\\}}}\Biggr{]}$
(56) $\displaystyle+$ $\displaystyle\frac{2\alpha
r^{2}\Bigr{[}-\frac{6M}{r^{4}}+\frac{4Q^{2}}{r^{5}}-\frac{m^{2}cc_{1}}{2r^{2}}+\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{r^{3}}\Bigr{]}^{2}}{\Bigr{[}1+4\alpha\Bigl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\Bigl{\\}}\Bigr{]}^{\frac{3}{2}}}$
$\displaystyle+$
$\displaystyle\frac{l(-2c^{2}c_{2}m^{2}r^{2}-m^{2}cc_{1}r^{3}+4Q^{2})}{r^{2}\sqrt{\Bigl{\\{}(-4c^{2}c_{2}m^{2}r^{2}-2m^{2}cc_{1}r^{3}+8Mr-4Q^{2})l^{2}-4r^{4}\Bigl{\\}}\alpha+r^{4}l^{2}}}.$

Therefore, using (55) and (56), Ricci and Kretschmann scalar are given by

$\displaystyle R$ $\displaystyle=$
$\displaystyle\frac{6}{\alpha}\Biggr{[}-1+{\sqrt{1+\frac{8\alpha
M}{r^{3}}-\frac{4\alpha Q^{2}}{r^{4}}-\frac{4\alpha}{l^{2}}-\frac{4\alpha
m^{2}c^{2}c_{2}}{r^{2}}-\frac{2\alpha m^{2}cc_{1}}{r}}}\Biggr{]}$ (57)
$\displaystyle-$
$\displaystyle\frac{\bigl{(}-4c^{2}c_{2}m^{2}r^{2}-m^{2}cc_{1}r^{3}+12Mr-8Q^{2}\bigl{)}^{2}\alpha
l^{3}}{2r^{2}\Bigr{[}\bigl{\\{}(-4c^{2}c_{2}m^{2}r^{2}-2m^{2}cc_{1}r^{3}+8Mr-4Q^{2})l^{2}-4r^{4}\bigl{\\}}\alpha+r^{4}l^{2}\Bigr{]}^{\frac{3}{2}}}$
$\displaystyle-$
$\displaystyle\frac{l\bigl{(}-10c^{2}c_{2}m^{2}r^{2}-3m^{2}cc_{1}r^{3}+24Mr-12Q^{2}\bigl{)}}{r^{2}\sqrt{\bigl{\\{}(-4c^{2}c_{2}m^{2}r^{2}-2m^{2}cc_{1}r^{3}+8Mr-4Q^{2})l^{2}-4r^{4}\bigl{\\}}\alpha+r^{4}l^{2}}},$
$\displaystyle K$ $\displaystyle=$
$\displaystyle\Biggr{[}\frac{1}{\alpha}\left[{1-\sqrt{1+4\alpha\left\\{\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\right\\}}}\right]$
(58) $\displaystyle-$
$\displaystyle\frac{4r\Bigr{[}-\frac{6M}{r^{4}}+\frac{4Q^{2}}{r^{5}}-\frac{m^{2}cc_{1}}{2r^{2}}+\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{r^{3}}\Bigr{]}}{\sqrt{1+4\alpha\biggl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\biggl{\\}}}}$
$\displaystyle+$ $\displaystyle\frac{2\alpha
r^{2}\left[-\frac{6M}{r^{4}}+\frac{4Q^{2}}{r^{5}}-\frac{m^{2}cc_{1}}{2r^{2}}+\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{r^{3}}\right]^{2}}{\Bigr{[}1+4\alpha\Bigl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\Bigl{\\}}\Bigr{]}^{\frac{3}{2}}}$
$\displaystyle-$
$\displaystyle\frac{r^{2}\Bigr{[}\frac{24M}{r^{5}}-\frac{20Q^{2}}{r^{6}}+\frac{2m^{2}cc_{1}}{r^{3}}-\frac{3m^{2}(2c^{2}c_{2}+cc_{1}r)}{r^{4}}\Bigr{]}}{\sqrt{1+4\alpha\Bigl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\Bigl{\\}}}}\Biggr{]}^{2}$
$\displaystyle+$
$\displaystyle\frac{4}{{r^{2}}}\Biggr{[}\frac{r}{{\alpha}}\left[1-\sqrt{1+4\alpha\left\\{\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\right\\}}\right]$
$\displaystyle-$
$\displaystyle\frac{r^{2}\Bigr{[}-\frac{6M}{r^{4}}+\frac{4Q^{2}}{r^{5}}-\frac{m^{2}cc_{1}}{2r^{2}}+\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{r^{3}}\Bigr{]}}{\sqrt{1+4\alpha\Bigl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\Bigl{\\}}}}\Biggr{]}^{2}$
$\displaystyle+$
$\displaystyle\frac{1}{{\alpha^{2}}}\Biggr{[}1-\sqrt{1+4\alpha\Bigl{\\{}\frac{2M}{r^{3}}-\frac{Q^{2}}{r^{4}}-\frac{1}{l^{2}}-\frac{m^{2}(2c^{2}c_{2}+cc_{1}r)}{2r^{2}}\Bigl{\\}}}\Biggr{]}^{2}$

(a) $\alpha=0.5$ (b) $\alpha=0.8$

Figure 20: $M=5$, $Q=1$, $c=1$, $c_{1}=-1$ and $c_{2}=1$. $m=0$ denoted by
solid black line, $m=1$ denoted by dash red line, $m=2$ denoted by blue dash
line, $m=3$ denoted by gold dash line,

In Fig. 20, Kretschmann scalar is plotted. The figure shows that $4D$
Einstein–Gauss–Bonnet massive gravity black hole has a true singularity at
$r=0$.