Charged Black Hole Solutions in Gauss-Bonnet-Massive Gravity

Motivated by high interest in the close relation between string theory and black hole solutions, in this paper, we take into account the Einstein-Gauss-Bonnet Lagrangian in the context of massive gravity. We examine the possibility of black hole in this regard, and discuss the types of horizons. Next, we calculate conserved and thermodynamic quantities and check the validity of the first law of thermodynamics. In addition, we investigate the stability of these black holes in context of canonical ensemble. We show that number, type and place of phase transition points may be significantly affected by different parameters. Next, by considering cosmological constant as thermodynamical pressure, we will extend phase space and calculate critical values. Then, we construct thermodynamical spacetime by considering mass as thermodynamical potential. We study geometrical thermodynamics of these black holes in context of heat capacity and extended phase space. We show that studying heat capacity, geometrical thermodynamics and critical behavior in extended phase space lead to consistent results. Finally, we will employ a new method for obtaining critical values and show that the results of this method are consistent with those of other methods.

may be modified to satisfy this point, where massive gravity is one such modification. From the theoretical perspective, the shear difficulty of constructing a consistent theory of massive gravity makes the subject more interesting. Therefore, there are limited considerable works in the massive gravity context. The development of the ghost-free theory with massive gravitons which are non-interactive in flat background was done in Ref. [20]. It was shown that generalization of this theory to curved background leads to existence of the ghost instabilities [21]. The effects of quantum interactions of massive gravity and a nonlinear theory of massive gravity in absence of ghost field [22] were investigated in Refs. [23,24]. In addition, a special class of charged massive black holes has been investigated [25] (see [26] for more details regarding considering massive gravity). A new class of nontrivial massive black holes in AdS spacetime was investigated in [27,28]. In this theory of the massive gravity, graviton has similar behavior as lattice in holographic conductor model. In other words, due to exhibiting a Drude peak which approaches to delta function in limit of massless gravity, graviton in this theory plays the role of lattice. Thermodynamical behavior, P − V criticality and geometrothermodynamics (GTD) of the Vegh's massive gravity has been, recently, investigated [29][30][31]. Some holographic consequences of graviton mass have been investigated in Ref. [32]. In this paper, we intend to generalize Einstein gravity action by adding GB Lagrangian and obtain exact charged solutions of GB-Massive gravity. In general, higher derivative gravity terms (such as GB) are corrections that become relevant at high-energy regime while the mass terms (of massive gravity) are relevant at low-energy. Here, we do not regard the energy regime, directly, and we look for the weight of GB and mass terms contributions into the thermodynamical behavior of a typical black hole.
Obtaining the solutions of GB-massive gravity, we should check their stability. In general, one may categorize the stability criteria into two classes; dynamical stability and thermodynamical one. In the present paper, we relax the dynamical stability and instead, we focus on thermodynamic behavior. Thermodynamical aspects of the black holes are among the popular subjects for studying black hole's behavior. Thermal stability of the black hole is of great importance and has been widely investigated in literature [33]. This is due to the fact that instability of black holes puts restriction (conditions) for validation of the physical solutions. One of the approaches for studying thermal stability is through canonical ensemble with investigation of the heat capacity. On the other hand, one can look for the phase transition points of solutions by calculating root(s) and divergence point(s) of the heat capacity. Also, the thermal stability and its corresponding conditions can be employed in studies conducted in context of AdS/CFT correspondence.
Recently, there has been a renewed interpretation regarding cosmological constant as a thermodynamical variable. It was pointed out that considering cosmological constant as a thermodynamical variable may lead to removing ensemble dependency in studying stability in context of canonical (heat capacity) and grand canonical (Hessian matrix) ensembles for BTZ black holes [34]. In addition, in context of AdS/CFT correspondence, it was shown that variation of cosmological constant in black holes corresponds to variation of the number of colors in Yang-Mills theory residing in the boundary of the spacetime [35]. There has been a consideration of the cosmological constant as a state-dependant parameter in two dimensional dilaton gravity [36]. In an interesting paper, Hawking and Page showed an interpretation of confinement/deconfinement phase transition in the dual strongly coupled gauge theory for phase transition between the stable large black hole and thermal gas in AdS space [37]. On the other hand, it was shown that there is a similarity between phase transition in black holes and liquid/gas phase transition in a Van der Waals system [38]. This consideration has been investigated in literature for different types of black holes [39].
Another method for studying critical behavior of the system is by constructing a spacetime by considering a thermodynamical potential and a set of extensive parameters as components of that thermodynamical spacetime. The calculated Ricci scalar of constructed thermodynamical spacetime should diverges in place of phase transition points that were mentioned in studying heat capacity. There are several approaches for this method that to name a few, one can say: Weinhold in which mass is considered as thermodynamical potential with other parameters (entropy, electric charge and etc.) as extensive parameters [40], Ruppeiner approach that is built up based on entropy as thermodynamical potential [41], it should be pointed out that it was shown that these two metrics are conformably related to each other by a conformal factor (temperature) [42]. The other method is the one that Quevedo introduced which has mass as thermodynamical potential like Weinhold but with different structure [43]. There are several types for Quevedo metric. Recently, it was shown that mentioned approaches in case of some black holes, will not yield consisting results with studying heat capacity. In order to overcome this problem, a new metric was introduced in Ref. [31,44]. The denominator of the Ricci scalar of HPEM metric [31,44] only contains expressions of numerator and denominator of the heat capacity.
The outline of the paper will be as follows. In Sec. II, we introduce action regarding charged massive gravity in presence of GB gravity and obtain charged black hole solutions in this gravity. We calculate conserved and thermodynamic quantities related to obtained solutions and check the validation of the first law of thermodynamics in Sec. III. In Sec. IV we study thermal stability of the massive charged GB black hole solutions in canonical ensemble using heat capacity. In the next section, considering cosmological constant as thermodynamical pressure, we study phase transition of black holes in context of P − V criticality and phase diagrams. In Sec. VI we study thermodynamical behavior of these solutions through geometrical thermodynamics. Sec. VII is devoted to calculation of critical values in extended phase space by using heat capacity. Last section will be closing remarks.

II. CHARGED BLACK HOLE SOLUTIONS IN GB-MASSIVE GRAVITY
The d-dimensional action of GB-massive gravity with a cosmological constant (Λ) can be written as where m and R are, respectively, the massive parameter and the scalar curvature and F = F µν F µν is the Maxwell invariant. In addition, α and L GB are, respectively, the coefficient and the Lagrangian of GB gravity, and Ψ is a fixed symmetric tensor. In Eq. (1), c i are constants and U i are symmetric polynomials of the eigenvalues of the d × d where R µν and R µνγδ are Ricci and Riemann tensors. Variation of the action (1) with respect to the metric tensor g µν and the Faraday tensor F µν , leads to where G µν is the Einstein tensor, H µν and χ µν are We consider the metric of d-dimensional metric with the following form where We should note that the constant k indicates that the boundary of t = constant and r = constant can be a positive (elliptic), zero (flat) or negative (hyperbolic) constant curvature hypersurface.
Here, we consider the ansatz of Ref. [29] for the auxiliary reference metric where c is a positive constant. Using the ansatz (7), U i 's may be written as [29] U 1 = d 2 c r , in which we used the notation d i = d − i. I order to study the effects of massive gravity, it is necessary to consider a reference metric, denoted by Ψ µν . In a special class, one may treat the reference metric in the context of a bimetric theory to credit it with a dynamical spin-2 tensor field (see [45] for more details). Another class of considering reference metric comes from the gauge invariance property with using of Stuckelberg scalar fields [46][47][48]. It was shown that the number of degrees of freedom propagated by any dRGT-type massive gravity in d-dimensional space-time (d ≥ 3) with N scalar fields (Stuckelberg fields) equals to 1 2 d(d − 3) + N . Using the Stuckelberg language it was shown that the number of propagators is less than the number of fields [49,50]. Vanishing of the determinant of the kinetic (Hessian) matrix of the scalar field Lagrangian helps us to reduce the number of propagators. Regarding the determinant of the kinetic matrix (detS), the so-called singular solutions are related to detS = 0 [51]. In general, we expect that all regular solutions (detS = 0) to be tangential to the surface detS = 0 at some point of time. It means that the singular solutions (detS = 0) are related to the families of regular solutions. In other words, we can set the conditions in vicinity of the singular surface and apply the infinitesimal evolution in time to obtain regular solutions which are tangential to the singular ones [52,53].
For example, 4-dimensional massive gravity theories admit, generally, up to six propagating modes where five of them are physical and one may be a ghost. Depending on the reference metric, the number of degrees of freedom that can be propagated, hence the existence of ghost may vary. The existence of the ghost determines whether the theory under consideration is actually stable or not. In this paper, the reference metric has the form of Ψ = (0, 0, ch ij ) where in a simpler case for 4-dimensional spacetime it will have a form of Ψ = (0, 0, 1, 1). This is a degenerate metric which is equivalent to a theory of massive gravity with Stuckelberg fields in the unitary gauge [27]. The reason for such consideration is the fact that under a specific coordinate transformation the spatial mass term breaks in spatial dimensions while the general covariance is preserved in radial and temporal coordinates. Since this choice of the reference metric leads to preservation of the Hamiltonian constraint, one of the degrees of freedoms is removed. In addition, the diffeomorphism is not broken in this specific theory of massive gravity. This leads to elimination of another degrees of freedom. Therefore, 4-degrees of freedom exist for the propagation with this reference metric. Now, considering that only two Stuckelberg fields exist in the diffeomorphism invariant formulation of this theory, two degrees of freedom are absent which leads to absence of Boulware-Deser ghost [27,54]. The existence of ghost for this theory was investigated in Ref. [22]. The mentioned reference metric is a singular one. It was shown that this theory with any singular metric is a ghost free theory [54]. In addition, a case of reduced massive gravity with two Stuckelberg fields was investigated in Ref. [52]. The extension of the theory to a negative cosmological constant and a Maxwell field was done in Ref. [53]. The stability of the black brane solutions and the degrees of freedoms that are provided by Maxwell field and massive gravity was investigated in detail and it was pointed out that the stability of the solutions and degrees of freedoms depends on the choices of the parameters. Regarding the above statements, we do not concern for the existence of regular solutions. Therefore in this paper we relax the dynamical stability criteria and focus on the thermodynamical properties as well as critical behavior. Use the electrical gauge potential ansatz, A µ = H(r)δ 0 µ , with Maxwell equation (3). we obtain where q is an integration constant which is related to the electric charge. Also, the Maxwell equation implies that the electric field in d-dimensions is given by Now, we want to obtain the static black hole solutions. For this purpose, one may use any components of Eq. (2) and obtain metric function f (r). One can use different components (tt and x 1 x 1 ) of Eq. (2) which can be written as where the prime and double prime are, respectively, the first and second derivatives with respect to r. We can obtain the metric function f (r) using the Eqs. (10) and (11), yielding where m 0 is an integration constant which is related to the total mass of the black hole. It is notable that, obtained metric function (12), satisfy the all components of the Eq. (2). Also, in the absence of massive parameter (m = 0), the solution (12) reduces to which describes a d-dimensional asymptotically (A)dS topological black hole with a negative, zero or positive constant curvature hypersurface. In order to study the effects of the massive gravity on metric function, we have plotted following diagrams (Figs. 1 and 2). It is evident that for specific values of different parameters, the metric function will have the usual behavior that one can see in GB-Maxwell gravity: there may be two horizons, one extreme horizon and no horizon (naked singularity) ( Fig. 1). Interestingly, considering specific set of values for metric function parameters will lead to modification in metric function behavior and variation from the usual behavior. In other words, the black holes in this configuration may have following behavior: two normal and one extreme (outer) horizons, four horizons (Fig. 2). The variation in number of the horizons emphasizes the contribution of the massive part. In other words, due to existence of the massive part, the geometrical structures of the black holes are modified and hence, their corresponding phenomenology is also altered. It should be pointed out that having three or four horizons is a function of the massive coefficients. Now, we are in a position to study the existence of the singularity. It is straightforward to show that all curvature invariants of the spacetime such as the Weyl square, Ricci square, Ricci and Kretschmann scalars are only the functions of f ′′ (r), f ′ (r) r and f (r) r 2 , thus it is sufficient to study one of them. It is easy to show that where in which Since the metric function diverges at r = 0 and is a smooth regular function for r > 0, its first and second derivatives are regular functions for r > 0. Therefore, we expect to have only one physical curvature singularity located at r = 0. To confirm it, one can investigate the Kretschmann scalar for various ranges of r. Straightforward calculations show that R αβγδ R αβγδ ∝ r −2(n−1) for small values of radial coordinate. In other words, one concludes which shows that there is an essential singularity located at the origin. For asymptotical behavior of the metric solutions, one can series expand the Kretschmann scalar for large values of radial coordinate which leads to in which for small values of GB parameter it will yield lim r−→∞ Eqs. (18) and (19) confirm that asymptotical behavior of the solutions is (A)dS with an effective cosmological Here we give a brief discussion regarding Carter-Penrose diagram. It is believed that for investigating the conformal structure of the solutions, one may use the conformal compactification method to plot the Carter-Penrose or conformal diagram (see Figs. 3,4 and 5). According to the Carter-Penrose diagrams one finds that, the singularity is timelike such as that of Reissner-Nordström black holes. In other words, although massive part of metric function may change the horizon structure of black holes, it does not affect the type of singularity and asymptotical behavior of the solutions. Drawing the Carter-Penrose diagrams, one can find that the causal structure of the solutions are asymptotically well behaved. For completeness, we should note that the curvature singularity comes from the nature of metric function and the geometrical nature of spacetime. In order to investigate physical singularity, one should study the Riemann tensor. It is a matter of calculation to show that the nonzero components of Riemann tensor for topological black holes are , and the first and second derivatives of the metric function are given before. Although we cannot analytically check the finiteness properties of the Riemann tensor, numerical calculations show that all singularities of the metric function (r = 0 and the roots of f (r)) are coordinate singularity except r = 0. The Kretschmann scalar and also the nature of Carter-Penrose diagrams confirm that the only physical singularity is located at r = 0.

III. THERMODYNAMICS
In this section, we calculate the conserved and thermodynamics quantities of the static black hole solutions in d-dimensional GB-massive context and then examine the first law of thermodynamics.
The Hawking temperature of the black hole at the outer horizon r + , may be obtained through the use of the definition of surface gravity. Straightforward calculations lead to T = 1 In order to calculate the electric charge of the black hole, one can use the flux of the electromagnetic field at infinity, yielding Next, for the electric potential, U , we use the following definition [12] In order to obtain the entropy of the black holes, due to generalization to Gauss-Bonnet gravity, one can use Wald's formula to calculate it [12]. This leads to which shows that area law is violated for GB black holes with non-flat horizons (k = 0). In order to obtain total mass of the black holes, one can use Hamiltonian approach which results into Having conserved and thermodynamic quantities at hand, we are in a position to check the first law of thermodynamics for our solutions. We obtain the mass as a function of the extensive quantities S and Q. One may then regard the parameters S and Q as a complete set of extensive parameters for the mass M (S, Q) where A is We define the intensive parameters conjugate to S and Q. These quantities are the temperature and the electric potential The results of Eq. (26) coincide with Eqs. (20) and (22) and, therefore, we find that these conserved and thermodynamic quantities satisfy the first law of black hole thermodynamics with the following form

IV. HEAT CAPACITY AND STABILITY IN CANONICAL ENSEMBLE
In this section, we study the stability of the solutions. In the context of black hole physics, there are two types of stability examination that one can employ, the so-called dynamical and thermodynamical stabilities. In this paper, we are only interested in thermodynamical stability of the solutions. We investigate thermal stability in context of canonical ensemble by calculating the heat capacity. Stability condition states that in order to a black hole being thermally stable, its heat capacity must be positive. There are two cases that may happen for an unstable black hole: it may have a phase transition and acquire stable state or the obtained solution is not a physical one (no phase transition takes place and the system is always unstable). Another important property of the heat capacity that motivates one to study it, is phase transition. It is stated that roots and divergence points of the heat capacity represent type one and type two phase transition, respectively. In other words, one can employ the numerator and denominator of the heat capacity for finding different phase transition points.
One can use following relation for calculating the heat capacity Considering Eqs. (20) and (23), it is a matter of calculation to show that heat capacity will be where As one can see, GB parameter and the topology factor (k) are coupled with each other. It shows that in case of flat horizon (k = 0), the effect of GB gravity is vanished. In other words, in horizon flat solutions, the heat capacity is independent of the GB gravity effects. In order to study the thermodynamical behavior of the system, stability and phase transition points, we have plotted Figs. 6 -10. It is evident that in case of spherical horizon (k = 1), there is a root for temperature, r R , in which for r + < r R , temperature is negative. Therefore, in this region, obtained solutions are non-physical. r R is a decreasing function of massive parameter (Fig. 6 right panel), an increasing function of electric charge (right panels of Figs. 9) and finally independent of GB parameter (Fig. 7 right panel). It should be pointed out that interestingly, increasing and decreasing mentioned parameters, will lead to formation of a maximum and minimum. Considering Eq. (28), these extrema are places in which heat capacity will have divergencies. In case of the topology, interestingly, for k = −1, temperature have a root and one divergence point (Fig. 8 right panel) which is not observed in case of horizon flat ( Fig. 8 middle panel) and spherical symmetric (Fig. 8 right panel).
For small values of massive parameter, the heat capacity is an increasing function of horizon radius ( Fig. 6 left panel). For specific values of the massive parameter, a maximum will be formed which is highly sensitive to variation of this parameter. For sufficiently large enough of m, there will be two divergencies for the heat capacity and between these two divergencies heat capacity is negative, hence, the system is unstable (Fig. 6 left panel).
Interestingly, the variation of the GB parameter has opposite effect on heat capacity comparing to variation of massive parameter. In other words, for small values of the GB parameter, black holes enjoy two divergencies with being unstable between these divergencies (Fig. 7 left panel). The interval between these divergencies is a decreasing function of GB parameter. For sufficiently small values of α, the heat capacity is only an increasing function of the horizon radius without divergency (Fig. 7 left panel).
As for topological effects, in case of spherical symmetric, no divergency is observed and for the r + > r R , system is in physical stable state (Fig. 8 left panel). On the other hand, in case of horizon flat, two phase transitions are observed for heat capacity with black holes being unstable between two divergencies (Fig. 8 left and middle panels). Finally, for k = −1, two roots and one divergence point are found. Between these roots and after divergence point, system is physical and stable whereas between larger root and divergence point, black holes are thermally unstable (Fig. 8 left and middle panels).
Similar behavior to variation of GB parameter is observed for variation of the charge except for one difference. For large values of GB parameter, heat capacity is only an increasing function of the horizon radius whereas by increasing the electric charge, the heat capacity will have maximum again ( Fig. 9 left panel).
Finally, in case of dimensions, it is observed that higher dimensionality has contribution to number and place of the divergencies for heat capacity. For 5-dimensions, only one root is observed, whereas in case of 6 and 7-dimensions, heat capacity has one root and two divergencies which are increasing functions of dimensions (Fig. 10). Now we are in a position to discuss phase transitions of the solutions. As it was stated before, roots and divergencies of the heat capacity are representing phase transition points of the system. In addition, it should be pointed out that thermodynamical principles state that systems in unstable states go under phase transition to stabilize. Therefore, for systems with one root and two divergence points, the following phase transitions happen: for the root there is a phase transition from unstable non-physical system to stable and physical one, for smaller divergence point a phase transition from larger black holes to smaller ones takes place and finally for larger divergency, systems go under phase transition from smaller black holes to larger ones.

V. P − V CRITICALITY OF CHARGED BLACK HOLES IN GB-MASSIVE GRAVITY
In this section, we study the phase transition points of the charged black holes in GB-Massive gravity through use of P − V criticality and related phase diagrams in spherical symmetric spacetime (k = 1). To do so, we consider following relationship between thermodynamical pressure and cosmological constant From thermodynamical point of view, one can indicate that conjugating thermodynamical variable corresponding to pressure would be thermodynamical volume. Therefore, in order to calculate the thermodynamical volume of the Considering cosmological constant as thermodynamical pressure leads to interpretation of mass not only as internal energy but also as Enthalpy of thermodynamical system. This interpretation leads to following relation for the Gibbs free energy of the system where by considering Eqs. (20), (23) and (24) in which By regarding this interpretation and using Eqs. (24), (30) and (31), one can obtain volume as which is in agreement with the volume of the spacetime with spherical boundary and radius r + . Because of the relation between volume and radius of the black hole, we use horizon radius (specific volume) for obtaining critical values. To calculate critical values, one can employ the inflection point properties As for thermodynamical pressure, by using Eqs. (20) and (30), one can obtain the following equation of state Now, we are in a position to calculate critical horizon radius. To do so, one should employ Eqs. (35) and (36) which lead to following relation 4d 2 q 2 6d 7/2 α ′ + d 5/2 r 2 + r 4 It is evident that obtaining critical horizon is not possible, analytically. Therefore, we use numerical method for obtaining critical horizon radius and corresponding critical temperature and pressure. The results of this numerical evaluation are presented in tables 1 and 2.  Considering obtained critical values for horizon radius, temperature and pressure, one can plot their corresponding phase diagrams (P − r + , T − r + and G − T diagrams). For the economical reasons, we plot phase diagrams for two set of critical values (Figs. 11 and 12). In order to elaborate the effect of variation of GB and massive parameters on critical behavior of the system, we plot two sets of comparing diagrams (Figs. 13 and 14).
The appearance of characteristic swallow tail in G − T diagrams shows that obtained values are critical ones in which phase transition takes place . Also, in P − r + diagrams the existence of region of phase transition and critical behavior of the system for different factors of the critical temperature are representing phase transition taking place . As for T − r + , the formation of subcritical isobars which divide the phase states to three different cases (before, in and after phase transition), indicates that calculated critical pressure is indeed the one that system goes under phase transition in .
It is evident that, critical horizon radius and universal ratio of Pcrc Tc are decreasing functions of massive parameter whereas critical temperature and pressure are increasing functions of it (see table 1 for more details). Interestingly, critical temperature and pressure are highly sensitive to variation of massive parameter. Increasing massive parameter for even small value leads to a significant change in the critical temperature. This shows that for black holes to have phase transition, it must be heated more significantly. In other words, for acquiring stable state, significant amount of energy is needed in order to have phase transition. It is worthwhile to mention that critical horizon also decreases highly but not with similar rate as temperature changes.
Interestingly, except for universal ratio of Pcrc Tc , the effects of variation of GB parameter on critical values, are completely opposite of massive parameter. In other words, critical temperature and pressure are decreasing functions of α whereas critical horizon is an increasing function of it (see table 2 for more details). It is worthwhile to mention that the effect of GB parameter is not as significant as the effect of massive parameter. Although both massive and GB gravities are extensions for Einstein gravity, their contributions are opposite. The value of GB parameter is representing the strength of the gravitational field. Taking this fact into account, one can conclude that strength of gravity puts restriction on (weakens) the contribution of the massive gravity.
Finally, it should be pointed out that subcritical isobars are increasing and decreasing functions of α and m, respectively. Subcritical isobars are representing the region in which phase transition takes place (Figs. 13 and 14 middle panels). Therefore, their increments/reductions will decrease/increase single state region of one phase which in our case it is whether small or large black holes. Interestingly, the gap between energy of two different phases is an increasing function of both GB and massive parameters. The only differences is that Gibbs energy level of two phases are increasing functions of massive parameter (Fig. 13

VI. GEOMETRICAL THERMODYNAMICS FOR HEAT CAPACITY AND P − V CRITICALITY
In this section we will study thermodynamical behavior of these black holes through geometrical thermodynamics (GTs). As it was pointed out earlier, there are several approaches for constructing thermodynamical spacetime. In this paper, we follow HPEM method for studying phase transition points of the heat capacity and later another metric for phase transition points of the extended phase space. The HPEM metric has the following form [31,44] where M S = ∂M/∂S, M SS = ∂ 2 M/∂S 2 and χ i (χ i = S) are extensive parameters. For economical reasons, we only plot corresponding figures for variation of α and m for studying the heat capacity (Figs. 6 and 7). By employing Eqs. (21), (23) and (25) and HPEM metric (Eq. (38)), one can plot following diagrams (Figs. 15 and 16).
Studying GTs diagrams shows that Ricci scalar of considered metric diverges in places of phase transition points. In other words, divergencies of the Ricci scalar coincide with divergencies and roots of the heat capacity (compare Figs. 15 and 16 with 6 and 7). Taking a closer look at the divergencies of the Ricci scalar, one can see that the behavior of the Ricci scalar around divergence points can be categorized into three groups: there is a change of sign around divergence point in which correspondingly there is a change of sign in heat capacity, divergency toward +∞ where there is a phase transition from larger black holes to smaller ones and finally divergency toward −∞ in which system goes under phase transition of smaller to larger black holes. This property of the Ricci scalar enables one to recognize the type of phase transition without studying heat capacity.
In Ref. [55] a new thermodynamical metric for studying critical behavior in case of P − V criticality was introduced. It was also shown that GTs, P − V criticality and heat capacity will lead to consistent results. In other words, number and places of phase transitions in these three pictures are uniform and they yield same results. There are two forms of GTs metric for studying critical behavior in extended phase space that are given as [55] In this paper, we only consider Case I metric for studying GTs in extended phase space. Considering Eqs. (24), (29), (30) and (39) with obtained critical pressure in table 1, we plot following diagram (Fig. 17).
In plotted graphs, a divergency is observed for Ricci scalar which is due to existence of the root for heat capacity (Fig. 17 left and down-right panels). For pressure smaller than critical pressure, two divergencies for heat capacity and Ricci scalar observed ( Fig. 17 up right panel). This behavior is similar to that was observed in Fig. 6 middle panel for pressures smaller than critical pressure. In other words, similar to T − r + diagrams, for pressures smaller than critical pressure system will have more than one phase transition point. Next, for pressure being critical pressure only one divergence point exists ( Fig. 17 down middle). The place of divergency is exactly where critical horizon radius was found. In other words, in case of considering critical pressure, the place of phase transition in extended phase space and divergencies of heat capacity and Ricci scalar coincide with each other (compare Fig. 17 down middle with Fig. 6 middle). Finally, for pressures larger than critical pressure no divergence point for heat capacity and Ricci scalar was observed which is consistent with result of T − r + diagram in Fig. 6.
According to the results, the GTs metric proposed for extended phase space is providing an effective method which its results are consistent with thermodynamical concepts. Also, this metric is able to point out phase transition points which one can obtain through studying T − r + and P − r + diagrams. More importantly, it is evident that heat capacity (its divergencies), GTs (its divergencies) and phase diagrams (P − r + , T − r + and G − T ) lead to consistent results.

VII. HEAT CAPACITY AND CRITICAL VALUES IN THE EXTENDED PHASE SPACE
The final section of this paper is devoted to calculation of critical values in extended phase space by using heat capacity. In last section, it was seen that critical values in which phase transition takes place in extended phase space present themselves as divergencies in heat capacity. Motivated by this result, in Ref. [55] a new method for obtaining critical values was introduced which is based on denominator of the heat capacity.
In this method, one consider denominator of the heat capacity and rewrite it with relation between thermodynamical pressure and cosmological constant. Next, one solve the denominator of the heat capacity with respect to pressure. Obtained relation for pressure is different from the one that is obtained by using temperature. If obtained relation for pressure has maximum(s), then in the place of that maximum (horizon radius and pressure), phase transition takes place. In other words, instead of mentioned method for obtaining critical values in section V , one can easily study the maximums of obtained relation for pressure through denominator of the heat capacity. The advantages of this method are pointed out in Ref. [55].
In order to obtain mentioned relation for pressure, one should use denominator of Eq. (29) with Eq. (30) which leads to in which −4q 2 2d 7/2 α ′ + d 5/2 r 2 + r 4 + + cd 2 2d 5 α ′2 + d 9 α ′ r 2 + + d 3 r 4 + r 2d2 + . Considering mentioned values in tables 1 and 2, we plot following diagrams for Eq. (40), (Fig. 18). It is evident from plotted graphs that the maximum of Eq. (40) is exactly located at the place in which black hole goes under phase transition (compare Fig. (18) with two tables) with same critical pressure. In other words, the maximums of Eq. (40) for different values of parameters, are representing critical horizon radii and pressures in which phase transitions take place. Therefore, by employing this method, one is able to calculate critical horizon radius and pressure at the same time without going through trouble of obtaining relation for calculating critical horizon radius and other relations. Another important property of this method is the consistency of pressure's behavior with thermodynamical concept. In other words, for pressures smaller than critical pressure, two values for pressure of phase transition are observed (corresponds to two divergencies for heat capacity in Fig. 17 and T − r + diagram in Fig. 6) and for pressures larger than critical pressure no pressure for phase transition point is observed.
In conclusion, one can see that this method enables us to obtain critical values and thermodynamical behavior of the system faster and its properties are consistent with thermodynamical concepts.

VIII. CLOSING REMARKS
In this paper, we have studied charged massive black holes in Gauss-Bonnet gravity. We obtained metric function for this gravity and showed that due to contribution of the massive part, the geometrical structure of black holes which includes number of the horizons and their places, have been modified. Therefore, the phenomenology differers completely from usual charged Gauss-Bonnet black hole. Plotted Carter-Penrose diagrams showed that although massive part of metric function can change the horizon structure of black holes, it does not affect the type of singularity and asymptotical behavior of the solutions. Next, we obtained conserved and thermodynamical quantities and showed that for obtained values the first law of thermodynamics is valid.
We also investigated thermal stability of the solutions in context of the canonical ensemble. We showed that variation of different parameters, affects the stability conditions of the black holes. Interestingly, it was found that the variation of massive and GB parameters has opposite effects on stability and phase transitions of the solutions. Also, it was shown that dimensionality has strong contribution to type and number of phase transition which resulted into modification of stability conditions. In other words, in specific dimensions, the black hole may only enjoy one type of the phase transition whereas in the other dimensions, it may have two types of phase transition in several points.
Next, by considering cosmological constant as thermodynamical pressure, we study phase transition of these black holes in the context of extended phase space. It was shown that thermodynamical volume is independent of GB and massive extensions and it only depends on topology of the solutions. Plotted phase diagrams showed that obtained critical values are the ones in which phase transition takes place. It was pointed out that interestingly, the effects of massive and GB parameters are opposite of each other. The critical temperature and pressure were highly sensitive to variation of massive parameter. In opposite, the GB parameter, hence strength of gravitational field put restriction on this rapid growth of temperature which will be a controlling factor for massive included systems.
In addition, through the use of GTs method, phase transition of these black holes was investigated in context of heat capacity and extended phase space. It was found that both employed geometrical metrics for studying thermodynamical behavior of the system yield consisting results with heat capacity and extended phase space. In other words, the divergencies of the Ricci scalar for these metrics coincided with divergencies of the heat capacity and thermodynamical behavior of the system in context of phase diagrams. Therefore, these three pictures yield consisting results.
Finally, a new method that was introduced in Ref. [55] employed for calculating critical pressure and horizon radius. It was shown that the results of this method (using maximum of obtained relation for pressure) and study conducted in extended phase space were in agreement and obtained critical pressure and horizon radius were the same.
It will be interesting to study the holographic aspects of obtained solutions in context of superconductors. Also, it will be worthwhile to study casual structure and casuality conditions of these systems. In addition, investigation of the causal structure of massive gravity in cosmological background by using of interesting Carter-Penrose diagrams is attractive. Moreover, following the approach of [52,53], one can study the stability and physical singularities of the solutions with massive degrees of freedom. We left these problems for the future works.