Criticality and extended phase space thermodynamics of AdS black holes in higher curvature massive gravity

Considering de Rham-Gabadadze-Tolley theory of massive gravity coupled with (ghost free) higher curvature terms arisen from the Lovelock Lagrangian, we obtain charged AdS black hole solutions in diverse dimensions. We compute thermodynamic quantities in the extended phase space by considering the variations of the negative cosmological constant, Lovelock coefficients ($\alpha_{i}$) and massive couplings ($c_{i}$), and prove that such variations is necessary for satisfying the extended first law of thermodynamics as well as associated Smarr formula. In addition, by performing a comprehensive thermal stability analysis for the topological black hole solutions, we show in what regions thermally stable phases exist. Calculations show the results are radically different from those in Einstein gravity. Furthermore, we investigate $P-V$ criticality of massive charged AdS black holes in higher dimensions, including the effect of higher curvature terms and massive parameter, and find that the critical behavior and phase transition can happen for non-compact black holes as well as spherically symmetric ones. The phase structure and critical behavior of topological AdS black holes are drastically restricted by the geometry of event horizon. In this regard, the universal ratio, i.e. $\frac{{{P_c}{v_c}}}{{{T_c}}}$, is a function of the event horizon topology. It is shown the phase structure of AdS black holes with non-compact (hyperbolic) horizon could give birth to three critical points corresponds to a reverse van der Waals behavior for phase transition which is accompanied with two distinct van der Waals phase transitions. For black holes with spherical horizon, the van der Waals, reentrant and analogue of solid/liquid/gas phase transitions are observed.


Introduction
Einstein's General Relativity (GR, also known as Einstein gravity) has been astonishingly regarded as the most successful description of gravitation and a well supported by numerous experiments since it was proposed [1][2][3] (see also these reviews [4,5]). Theoretically, inconsistency appears when GR is supposed to be reconciled with the laws of quantum a e-mail: hendi@shirazu.ac.ir b e-mail: ali.dehghani.phys@gmail.com physics for producing the theory of quantum gravity. From the experimental side, Einstein gravity has a problem with the accelerated expansion of the universe in the large scale structure since it needs an unknown source of energy (the so-called dark energy) captured by the cosmological constant [6][7][8][9]. In this regard, various attempts have been made to find an alternative such that it modifies Einstein gravity in the large scales (IR limit). Massive gravity is one of the alternatives that modifies Einstein gravity by giving the graviton a mass and provides a possible explanation for the accelerated expansion of the universe without the requirement of dark energy component [10][11][12]. Assuming that gravitons are dispersed in vacuum like massive particles, gravitational waves' observation of the coalescence for a pair of stellarmass black holes (GW170104) has bounded the graviton mass to m g ≤ 7.7 × 10 −23 eV /c 2 [13]. On the other hand, depending on the exact model of massive gravity, the graviton mass is typically bounded to be a few times the Hubble parameter today, i.e., m g ≤ 10 −30 − 10 −33 eV /c 2 , in which for graviton mass region m g 10 −33 eV /c 2 , its observable effects would be undetectable [14] (for more details on different mass bounds see [15]). Massive gravitons, if they exist, are yet to be found; but, according to the recent data of LIGO, such an assumption is experimentally logical and therefore deserves to be explored theoretically [10][11][12][13][14][15].
Depending on what features of GR is accepted unchanged, various theories of gravity have been created. Modification of GR is characterized by a deformation parameter such as Lovelock coefficients α i 's in Lovelock gravity (which determines the strength of higher curvature terms) or graviton mass parameter in massive gravity models. Based on the nature of the deformation parameter, the original theory (GR) can be recovered by taking some limits (e.g. the zero limit of graviton mass parameter must recover GR and its associated outcomes). In order to have a generalized and well-defined theory, we should take care of ghosts. Although the first linear version of the massive theory (i.e., the Fierz-Pauli model [16]) is ghost-free, it does not lead to the linearized GR as the graviton mass goes to zero which is known as the vDVZ discontinuity [17,18]. Vainshtein discovered that such discontinuity appears as a consequence of working with the linearized theory of GR [19], and by employing the Stueckelberg trick it can be found that all degrees of freedom introduced by the graviton mass do not decouple in the zero limit of graviton mass [20]. On the other hand, Boulware and Deser showed some specific nonlinear models of massive gravity suffer from ghost instabilities however they could restore continuity with GR [21]. Eventually, the de Rham-Gabadadze-Tolley (dRGT) theory of fully nonlinear massive gravity resolved the ghost problem in four dimensions by adding higher order graviton self-interactions with appropriately tuned coefficients [22][23][24][25]. The higher dimensional extension of the massive (bi)gravity has been discussed in [26,27], which confirms the absence of ghost fields using the Cayley-Hamilton theorem. Interestingly, in massive gravity framework, spherically symmetric black hole solutions were found in [28,29] and in the limit of vanishing graviton mass they go smoothly to the Schwarzschild and Reissner-Nordström (RN) black holes. Furthermore, asymptotically flat black hole solutions were found in [30], but the curvature diverges near the horizon of these solutions. In this regard, black hole solutions with non-singular horizon were introduced in [31] with the identification of the unitary gauge to the coordinate system in which black hole has no horizon (for more details see [31]). The other interesting solutions related to the cosmology, gravitational waves and (time dependent) black holes were found in [10,[32][33][34][35][36][37][38] which will not discuss in this paper. Of interesting case for us is Vegh's black hole solution [39] in which the general covariance preserves in the "t" and "r " coordinates, but, is broken in the other spatial dimensions. In [40], this solution was generalized to topological black holes in higher dimensions. Inspired by the interesting features of these solutions, black hole solutions of massive gravity coupled to the higher curvature terms, dilaton and nonlinear electromagnetic fields were constructed and studied in detail [41][42][43][44][45][46][47].
From the string theory point of view and also braneworld cosmology perspective, Lovelock gravity [48,49], as a natural deformation of GR in higher dimensions [50], has an essential role. The most important motivation to study such a theory is related to superstring theory models which lead to ghost-free nontrivial gravitational interactions in higher dimensions [51]. The low-energy limits of type II string theory and E 8 × E 8 heterotic superstring give rise to effective models of gravity in higher dimensions which contain higher powers of the Riemann curvature (e.g., R 2 , R 3 , R μν R μν , R μνγ δ R μνγ δ , . . .) in addition to the usual Einstein and cosmological constant terms [52][53][54][55]. It is notable that the ghost-free combinations of these terms are proportional to the Euler invariant [51,56] which is exactly the same as the Lovelock Lagrangian [48,49]. The Lagrangian of the Lovelock gravity is given by a sum of dimensionally extended Euler densities as in which δ μ 1 ν 1 ... μ k ν k ρ 1 σ 1 ... ρ k σ k and R μ k ν k ρ k σ k are the generalized totally antisymmetric Kronecker delta and the Riemann tensor respectively. In 4-dimensional spacetimes, the Lovelock theory of gravity reduces to GR, and in spacetime dimensions with d 5 the Gauss-Bonnet term appears, and for d 7 the third order Lovelock term has contribution besides the Einstein and Gauss-Bonnet terms. Lovelock gravity is also ghost free with second order field equations and admits black hole solutions and the associated thermodynamics as expected [57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75].
Of interesting case for theoretical physicists is thermodynamic properties of black holes in comparison with ordinary systems in nature. In fact, black hole mechanics obeys the same laws as the laws of thermodynamics [76], and many investigations have confirmed this statement for more complicated black hole spacetimes in modified gravities. In addition, a wealth of results in the context of black hole thermodynamics have been presented which show black holes in Einstein gravity can imitate some thermodynamic properties of ordinary systems such as the van der Waals (vdW) phase transition which represents a liquid/gas (first order) phase transition [77], the reentrant phase transition (RPT) in multicomponent fluid systems [78,79] and the triple point in solid/liquid/gas phase transition. The phase space of Schwarzschild-anti de Sitter (Schwarzschild-AdS) black holes admits the so-called Hawking-Page phase transition [80] which is interpreted as a confinement/deconfinement transition in the dual boundary gauge theory (SYM plasma) [81]. Remarkably, RN-AdS and Kerr-Newman-AdS black holes possess a first order phase transition which closely resembles the well-known vdW phase transition in fluids [82][83][84]. Interestingly, Born-Infeld-AdS black holes, as a nonlinear version of RN-AdS ones, display a phase structure which relates the mass (M) and the charge (Q) of the black holes similar to the solid-liquid-gas phase diagram [85]. These considerations were done in the presence of the cosmological constant as a fixed parameter and recently are referred in the literature as non-extended phase space. In fact, as stated in [86], these mathematical analogies are confusing since some black hole intensive (extensive) quantities have to be identified with an irrelevant extensive (intensive) quantities in the fluid system, for example, the identification between the fluid temperature and the charge of the black hole is puzzling.
In thermodynamic systems, some quantities are thermodynamic variables and the others are fixed parameters which cannot vary. Only experiment can determine that a quantity (parameter) can vary or hold fixed. From the theoretical perspective, one can assume the variation of a fixed parameter of a theory and then see its consequences. In this regard, the later mismatch between extensive and intensive quantities of the black hole and fluid systems can be solved if one treats the cosmological constant ( ) as a thermodynamic variable, i.e., pressure [86]. This idea (which first established in [87] and then developed in [88][89][90][91][92]) leads to an extension of the phase space thermodynamics and the exact analogy between quantities of a black hole and liquid-gas systems at the critical point. For example, a transition occurs in P − T plane for the both of RN-AdS (small/large) black hole and liquid-gas systems. In addition, the variation of in the first law of black hole thermodynamics solves the inconsistency between the Smarr formula and the traditional form of first law since in the presence of a fixed cosmological constant the scaling argument [87] is no longer valid. This motivates consideration of the first law of black hole thermodynamics with varying which is referred to as the extended phase space thermodynamics in physics community. Regarding the extended phase space thermodynamics, RPT has been observed for Born-Infeld-AdS and singly spinning Kerr-AdS black holes in the context of Einstein gravity [93,94]. In a black hole system, it is interpreted as large/small/large black hole phase transition. Moreover, the analogue of solid/liquid/gas phase transition was found for doubly spinning Kerr-AdS black holes which is interpreted as small/intermediate/large black hole transition [95,96].
The objective of this paper is to construct the higher curvature massive gravity in order to study the effects of higher order curvatures on the black hole solutions of massive gravity and investigate the associated criticality and thermodynamics in the extended phase space. Indeed, some thermodynamic features of black holes, e.g. universal ratio, may depend on the specific choice of the gravitational theory. Therefore it is so important to understand the effect of modified gravities. We select the Lovelock gravity up to third order (referred to as TOL gravity) as the higher curvature framework for our investigations. When the Lovelock massive theory of gravity (LM gravity) is constructed, in principle, the parameters α (Lovelock coefficient) and m (graviton mass) are considered as deformations of GR, and by taking the limits m → 0 and α → 0, GR is naturally recovered. According to scaling argument, any dimensionful parameter in a given theory has a thermodynamic interpretation and as a result, the Smarr formula and the first law of black hole thermodynamics must be modified. According to this fact, thermodynamically, more interesting phenomena can take place in a more complicated theory of gravity such as Lovelock and massive gravities which have a finite number of dimen-sionful parameters. One can observe that modified gravities such as massive and Lovelock theories exhibit a rich black hole phase space structure with respect to those counterparts in Einstein gravity. The existence of higher order curvatures based on the third order Lovelock (TOL) gravity can lead to critical behavior and phase transition for AdS black holes with hyperbolic horizon topology [97][98][99][100] in contrast to Einstein gravity which only spherically symmetric AdS black holes admits phase transitions. Remarkably, hyperbolic vacuum black holes in Lovelock gravity expose non-standard critical exponents at a special isolated critical point which are different from those of vdW ones [101]. Until writing this paper, a wealth of evidence has been indicating that all the black hole solutions in Einstein gravity in the presence of any matter field have the same critical exponents as the vdW fluid [86,93,102,103]. Interestingly, a "λ-line" phase transition occurs for a class of AdS-hairy black holes with hyperbolic horizon in Lovelock gravity where a real scalar field is conformally coupled to gravity [104]. In addition, for charged black branes, the inclusion of higher curvature gravities based on a generalized quasi-topological class could lead to phase transition and critical behavior with the standard critical exponents [105]. These indications reveal the rich phase space structure of Lovelock gravity's black holes. On the other hand, in the massive gravity framework, phase transition and critical behavior could take place for all kinds of topological black holes [106]. In this regard, the vdW and RPTs were found for AdS black holes [107,108], and in the presence of Born-Infeld (BI) nonlinear electromagnetic fields, the triple point emerges and the corresponding large/intermediate/small transition could take place [109].
Taking these considerations seriously, in this paper, we mainly focus on the critical behavior and phase transitions of AdS black hole solutions in the Lovelock massive (LM) gravity. By constructing this model, besides its novel phase structure, we could be able to figure out what characteristic features of Lovelock and massive gravities persist or ruin. Thus, we have organized this paper as follows: first, in Sect. 2.1 , we give a brief review of thermodynamics in extended phase space, stability analysis and phase transition for AdS black holes in the context of Einstein gravity. After, in Sect. 3, we construct the LM gravity by introducing the action and associated field equations and then present a new class of charged-AdS black hole solutions in arbitrary dimensions. By computing thermodynamic quantities, we prove the traditional first law of black hole thermodynamics is satisfied. After that, in Sect. 4.1, we perform a thermal stability analysis for the obtained black hole solutions in the canonical ensemble. Furthermore, we reconsider the first law of thermodynamics in the extended phase space and then study P − V criticality and phase transition(s) of black holes to complete our discussion. Finally, in Sect. 5, we finish our paper with some concluding remarks.

General formalism: thermodynamics, stability and phase transition for AdS black holes
In this section, we will develop the basic framework that we need to study the critical behavior and thermodynamic properties of AdS black holes in the next sections. Some useful issues will be briefly reviewed such as gravitational partition function, black hole thermodynamics, local thermodynamic stability, phase transition and critical behavior of black holes. Throughout this paper, we use the geometric units, G N =h = c = k B = 1. In these units,

Basic set up: partition function and action
According to an old idea of unification, it is believed that all known forces (strong, weak, electromagnetic and gravitation) in nature might be unified in the so-called "theory of everything". For many years, physicists have been looking for a consistent theory in which all forces in nature are described using path integral formalism of quantum field theory (QFT), like QCD and electroweak theories. Respecting such approach, we expect generating functional of the quantum theory of gravity could be defined by a Euclidean path integral over a dynamical metric (tensor field), g μν , as follows where φ is considered as matter fields and I G represents the on-shell gravitational action which is obtained by substituting the classical solutions of g. The generating functional Z contains a complete summary of the theory which its dominant contribution originates from the classical solution of the action by applying the stationary phase method (also known as steepest descent method or saddle-point approximation). Since the Euclidean formalism is obtained by applying the Wick rotation (t E = it) on the Lorentzian version, the Euclidean spacetime would be periodic in time. Following the method proposed by Matsubara [110], one can use the mapping to calculate the partition function of a thermodynamic system applying the techniques of calculus in QFT, and vice versa (for more specific details see [111][112][113]). As a result, considering the substitution (2.2), it is natural to regard Z as the statistical mechanical partition function of a gravitational system such as a black hole. Comparing Eq. (2.1) with the free energy (F = −k B T ln Z), we find our main relation as It should be noted that F can be identified with Helmholtz or Gibbs free energy depending on the thermodynamic variables of the system. Therefore, thermodynamic quantities associated with the gravitational system can be directly extracted by using the standard methods in statistical mechanics [114,115]. For systems in nature, one can define an appropriate action in which its dynamical equation of motion is determined by the variational principle. The total action (I G ) for a gravitational system, here black holes, consists of three terms as where I b , I s and I ct are, respectively, called the bulk action, the surface term (boundary action), and the counterterm action. The surface term is needed to have well-defined variational principle and remove the derivative terms of g μν normal to the boundary, and the counterterm actually is an alternative surface term that comes from renormalization method in QFT to eliminate UV divergences and only can regulate asymptotically AdS spacetimes [116][117][118][119][120][121]. For asymptotically flat spacetimes, one can use the subtraction method to cancel divergences [122,123]. A finite number of surface terms and counterterms are always needed to have a set of well-defined field equations and a finite total action. Since we intend to study the critical behavior of black holes and understand the theory-dependence behind this phenomenon, we devote the rest of this paper to explore in these objects. Of interesting case is d-dimensional topological AdS black holes with the following metric ansatz where the line element h i j dx i dx j is the metric of (d − 2)dimensional (unit) hypersurface with the constant curvature d 1 d 2 k and volume ω d 2 with the following form in what follows we will use this notation). The different values of the topological factor (i.e., k = −1, 0, +1) determine the topology of the event horizon and could be positive (spherical, S n ), zero (Planar, R n ), or negative (hyperbolic, H n ). The details of metric function ψ(r ) depend on the theory that we pick out. In this paper, we always consider our line element ansatz as above. Now, we focus on the charged (static) AdS black holes in Einstein's GR, briefly. The bulk action for the Einstein gravity on the d-dimensional background (manifold M) in the presence of negative cosmological constant and (Maxwell) electromagnetic filed is where g is the determinant of metric tensor g μν , = −(d−1)(d−2)/2 2 with the AdS radius , and F = F μν F μν is the Maxwell invariant in which F μν = ∂ μ A ν −∂ ν A μ is the electromagnetic field tensor (Faraday tensor) with the gauge potential A μ . The Einstein bulk action has to be accompanied by boundary action(s) and counterterm(s). The Gibbons-Hawking boundary action and the counterterm for regulating divergences of the Einstein bulk action have been introduced in [124] and [119], respectively. In this paper, we leave out the details of these terms and refer the readers to the above references where the relevant details can be found.
Gravitational field equations are obtained by varying Eq. (2.7) with respect to the metric tensor g μν as where G μν = R μν − 1 2 g μν R is the Einstein tensor. Considering the line element ansatz (2.5), the so-called RN-AdS black hole solutions of gravitational field equations (2.8) are given by where q and m 0 are integration constants related to the electric charge and finite mass of the black holes, respectively. Calculating the Kretschmann scalar, one can find an essential singularity. In fact, the Kretschmann scalar diverges only at r = 0 which is covered with an event horizon and thus one can interpret it as a black hole. The larger root of ψ(r + ) = 0 with a positive slope determines the event horizon of black holes. In Einstein gravity, numerical calculations show that the metric function may have two real positive roots (RN black holes), one extreme root (extreme black holes) or it may be positive definite (naked singularity). These results still hold for higher curvature gravities such as the Lovelock gravity [70][71][72]74], but the situation is different in massive gravity for which the metric function may have more than three roots [41,44,45] as will be shown in Sect. 3.2.

Black hole thermodynamics
One can regard charged (static) black holes as thermodynamic systems with the following first law [76] d M = T d S + d Q, (2.10) in which a formal equivalence between the physical temperature T and the surface gravity (κ) was established in [125,126] by studying quantum fields and particle creation near the event horizon for different observers in the black hole background. The notion of black hole, like ordinary thermodynamical systems, should be assigned an entropy and a temperature. Entropy of the black holes satisfies the so-called area law (S = A/4) in the Einstein gravity and this may be modified for higher derivative gravities (see Refs. [70][71][72] for the modified entropy in the Lovelock gravity). In addition, respecting the first law of thermodynamics, for all gravitational theories, one can always obtain the entropy using the following relation which needs the functional forms of mass (M) and temperature (T ). Besides, one can use the definition of the surface gravity with the Killing vector χ as (2.12) to calculate the Hawking temperature as It is worth mentioning that χ = ∂ t is the temporal Killing vector for static spacetimes. In order to define the electric potential, it is necessary to select a reference. Naturally, electric potential can be measured at the horizon with respect to infinity as a reference, i.e., (2.14) Moreover, the electric charge is an extensive quantity corresponds to the potential as an intensive quantity. Using the Gauss' law as and calculating the flux of the electromagnetic field at infinity, the electric charge is obtained. Finally, the mass M of the static black hole can be written down as measured by a faraway observer using Ashtekar-Magnon-Das (AMD) formula as [127][128][129] in which m 0 is a positive integration constant (see Eq. 2.9) and easily obtained from ψ(r + ) = 0 (its details depend on the parameters of the gravitational theory). This formula still holds for higher curvature gravities like Lovelock theory [70,72,74] and can be calculated using the behavior of the metric at large r (for asymptotically flat black holes) or counterterm method (for asymptotically AdS black holes).
In the following Obviously, in the traditional treatment of the first law of black hole thermodynamics, the work term "PdV " is missed. Recent developments [86][87][88][89][90][91] show that one can extend the thermodynamic phase space and insert the volume-pressure term in the first law by treating the negative cosmological constant ( ) as a thermodynamic variable. If one employs the following identification naturally, its conjugate variable will be interpreted as the thermodynamic volume. It should be noted that, here, the gravitational background necessarily is an asymptotically AdS spacetime in order to have a positive definite pressure. Such identification leads to an extended version of the first law of black hole thermodynamics with a volume-pressure term. In the extended phase space, the first law is obtained as Since the above first law has been written in the enthalpy representation (the Legendre transform of the energy representation, H = E + PV ), the mass of black hole has to be interpreted as the enthalpy (M ≡ H ). Consequently, the thermodynamic volume is calculated as in which X i 's are the extensive quantities Q and S. Now, it can be tested if the thermodynamic quantities of RN-AdS black holes satisfy the first law. Thermodynamic quantities for RN-AdS black holes read

Thermal stability of black holes
Here, we are going to explain the basic treatment for studying the thermal stability of black holes. Analyzing the local thermal stability of a black hole can be performed in both canonical and grand canonical ensembles. It can be done by studying the behavior of entropy S near a sufficiently small neighborhood of a given point in the phase space of possible extensive variable X i . Entropy is a smooth function of extensive variables in which they are M and Q in our case. Local thermal stability states that entropy must be a smooth concave function of the extensive variables (S = S(M, Q)). This is equivalent to have a negative definite value for the determinant of the Hessian matrix of the entropy S, i.e. [130]. Also, since one could invert S = S(M, Q) picture to M = M(S, Q) one, the stability criterion H S X i ,X j ≤ 0 may be expressed, differently, as the determinant of the Hessian matrix of M being positive definite [131], i.e.
In what follows, we will use this later criterion in order to analyze thermal stability of AdS black holes.
In the canonical ensemble, the total charges are held fixed and, consequently, the Hessian matrix has one component as H M S,S ≡ ∂ 2 M/∂ S 2 . As we have shown below, the Hessian matrix will be the function of the heat capacity C Q . It is possible that we confront with a situation for which both of T and H M S,S are negative definite, and therefore, C Q acquires a positive value. In order to avoid this problem one always has to take care of positivity of the temperature and heat capacity, simultaneously. In conclusion, positivity of the heat capacity ensures the local thermal stability in the nonextended phase space of allowed physical black hole quantities with T, M 0.
For extended phase space thermodynamics in the canonical ensemble, the pressure as well as total charges is kept constant. Consequently, the stability criterion (2.26) naturally leads to positivity of heat capacity at constant pressure and electric charge, i.e., where M ≡ H . Of course, one can still define another ensembles in which pressure or U (1) charge could vary. In our considerations, we will perform stability analysis in the extended phase space of the canonical ensemble. In this regard, the results are qualitatively the same as the case of non-extended phase space, but thermodynamic interpretations are basically different. Now, we perform a thermal stability analysis for RN-AdS black holes. To determine thermally stable regions for black holes, first, one has to find in which regions the associated temperature is positive. That could depend on the topology of event horizons. In Fig. 1, typical behaviors of temperature are depicted for different horizon topologies (k = 1, 0, −1) in four and higher dimensions. Clearly, there is always a bound point for the radius of the event horizon (r b ) in which the temperature of the black hole is positive for r + > r b . The topological type of event horizon (k = 1, 0, −1) can change the value of r b . The more value for k results into the lower value for r b . In addition, in the region r + > r b the mass of the black holes (M) is always positive.
The heat capacity of RN-AdS black holes can be calculated as Using the heat capacity, one can obtain some information about phase transition. For example, a divergence point of heat capacity is a sign of a possible phase transition. For large enough values of r + , since temperature and heat capacity are both increasing functions of r + , thus the large black holes are thermally stable. Here, we examine the heat capacity for RN-AdS black holes with different topological factor (k) in more details. This will complete our discussion of local thermal stability.
Spherical black holes (k = 1): In this case, the heat capacity has only one positive root (r b ) in which C Q is negative definite for regions r + < r b . For regions r + > r b , depending on the values of q and spacetime dimensions (d), three possibilities may happen: (i) C Q is an increasing function of r + , thus, in regions r + > r b , the heat capacity will be positive and black holes are thermally stable. (ii) C Q may have two divergence points for the RN-AdS black holes, where between divergence points the heat capacity is negative definite (unstable black hole region), thus a thermal phase transition can happen. (iii) C Q has one divergence point which is positive around such divergency. Such a single divergence point, with positive C Q around it, may indicate critical behavior of the system. In Fig. 2, the possibilities of items (i) and (ii) are depicted. We refer to the first and second divergence points as r m and r u respectively (r b < r m < r u ). In regions r b < r + < r m and r + > r u , RN-AdS black holes are thermally stable (a phase transition could occur between these two thermally stable regions) and for the other regions are unstable.
Ricci flat black holes (k = 0): In this case, the heat capacity has only one positive root, r b , in which C Q is negative definite for regions r + < r b and positive for r + > r b . According to Eq. 2.28, the heat capacity does not diverge for finite values of r + since the denominator of the heat capacity cannot have any root (see the left panel of Fig. 3). As a result, Ricci flat black holes are thermally stable for regions r + > r b and no phase transition takes place. Hyperbolic black holes (k = −1): In this case, as well as Ricci flat black holes (k = 0), there is only one root (r b ) for the heat capacity. Again, hyperbolic black holes are thermally stable in regions r + > r b since the heat capacity is always positive and no phase transition takes place (see the middle and right panels of Fig. 3).

Phase transition for RN-AdS black holes
In this section, we review the essence of critical behavior and phase transition in AdS black holes. In recent years, the interesting analogy between liquid/gas and small/large black hole phase transitions has attracted the attention of many authors. In this regard, the exact analogy between liquid-gas system (vdW fluid) and charged-AdS black holes was first completed by [86] in the context of Einstein-Maxwell gravity. In fact, RN-AdS black holes exhibit first-order phase transition with the same critical exponents as the vdW system. We will gen-eralize the results in [86] for higher dimensions and various event horizon topologies and will show that P − V criticality only exists for spherically symmetric black holes.
Here, we show the extension of thermodynamic phase space, by introducing negative cosmological constant as thermodynamic pressure (i.e. P = − /8π ), indicates the phase transition for charged-AdS black holes. Using Eq. (2.20), the equation of state reads as in which P phys and T phys denote the physical pressure and temperature, respectively. As a result, physical equation of state is given as Equation (2.30) should be compared with the following equation of state of vdW fluid in which a (molecular interaction forces) and b (molecular size) are positive definite quantities, and v is the specific volume of vdW fluid. To do so, rewriting vdW equation in terms of P and then expanding for v b lead to Comparing Eq. (2.31) with the above expansion implies that the horizon radius (not the thermodynamic volume V ) is associated with the vdW fluid specific volume, i.e.
It should be emphasized, in this analogy, the same physical quantities are compared with each other. To sum up, we summarize the analogy between vdW fluid and RN-AdS black hole in the following table.
vdW fluid RN-AdS black hole The critical point occurs at the spike like divergence of specific heat at constant pressure (i.e., an inflection point in the P −V (P −r + ) diagram) and can be found by considering the following equations, simultaneously The thermodynamic quantities P c , T c and v c must be positive definite. Evidently, for uncharged-AdS black holes (q = 0) or non-spherical horizon solutions (k = 0, −1), there are no phase transition and critical behavior. In fact, the electric charge parameter q and topological factor k play crucial rules to have critical behavior for AdS black holes in the Einstein gravity. In Sect. 4.1, we will show this is not the case in Lovelock or massive gravities and phase transition can happen even for uncharged black holes or black holes with Ricci flat/hyperbolic horizons. The thermodynamic quantities at the critical point satisfy the following universal ratio and of course, is only valid for spherical black holes (k = 1). Amazingly, the universal ratio for RN-AdS black holes coincides with the vdW fluid (the universality is equal to 3/8 in 4-dimensions) and it only depends on spacetime dimensions. Now, by finding and analyzing the free energy of RN-AdS black hole system, we complete our discussion on the criticality. In the fixed charged ensemble (canonical ensemble), the on-shell action is identified with the Gibbs free energy (since is a thermodynamic variable in the theory). When Q is held fixed, one has to add a surface term (for electromagnetic field) to fix charge on the boundary. As a result, the total action (2.4) has to be accompanied with the following boundary term The Gibbs free energy can be obtained using the Legendre transformation or calculating the on-shell action as follows (2.39) The qualitative behavior of Gibbs free energy as a function of temperature is depicted in the right panel of Fig. 4. Obviously, the swallow-tail behavior demonstrates the first order phase transition (exactly the same as vdW fluid). For the sake of completeness, in the left and middle panels of Fig. 4, P − r + and T − r + diagrams are plotted. Evidently, the RN-AdS equation of state (2.29) mimic the behavior of the vdW fluid for any fixed Q.

Action and field equations
In this section, we start the investigation of black holes in the LM gravity framework. Again, the total action consists of Middle and right panels: P < P c (continuous lines), P = P c (dotted lines) and P > P c (dashed lines) three terms, I G = I b +I s +I ct . The Lovelock boundary term (I s ) is needed to have a well-defined gravitational action, and the counterterm for static solutions with Ricci flat or curved horizons (I ct ) is necessary to regulate divergences of the Lovelock bulk action. The boundary and counterterm of Lovelock gravity are constructed before and the relevant details can be found in [71,132,133]. We add generic mass terms for the gravitons in the Lovelock theory of gravity by supplementing higher order interaction terms of massive gravitons (which first proposed in [23]). Considering this theory on the d-dimensional background manifold M with the dynamical metric g μν , the bulk action for LM gravity with a U (1) gauge field (Maxwell field) may be written as where α 0 L 0 = −2 and α 1 L 1 = R are, respectively, the negative cosmological constant and the Ricci scalar, and L 2 and L 3 are called the Gauss-Bonnet (GB) and the third order Lovelock (TOL) Lagrangians in the literature which have the following forms The Lovelock coefficients α 2 and α 3 , which are positive definite constants [57], indicate the strength of the second and third order curvature terms. In the action (3.1), the last term is called massive interaction terms in which m is the graviton mass parameter, and f is the auxiliary reference metric which is a fixed rank-2 symmetric tensor. The massive coupling c i 's, as the free parameters of the theory, are arbitrary constants which their values can be determined according to observational or theoretical considerations [10][11][12]31,39].
g μα f αν with the following explicit forms [22,23,39] where the square root in K stands for matrix square root, , and the rectangular bracket denotes the trace [K] = K μ μ . We restrict our study to U i up to the fourth interaction term (U 4 ), while considering the higher order terms is straightforward.
Using the variational principle, the electromagnetic and gravitational field equations of LM-Maxwell gravity can be obtained as in which G μν = R μν − 1 2 g μν R is the Einstein tensor and G (GB) μν and G (TOL) μν are, respectively, the second (GB) and third order Lovelock (TOL) tensors given as In addition, X μν is given by In the next section, we will obtain charged-AdS black hole solutions of the fully nonlinear gravitational field equations (3.4) and study their geometric properties.

LM charged-AdS black holes
Here, we intend to obtain static black hole solutions of LM gravity. In order to achieve the topological charged-AdS black holes, we consider the d-dimensional line element ansatz given in Eq. (2.5). We make use of the appropriate ansatz for the reference metric f μν with the following form [39,40,45] where c 0 is a positive constant. It is important to note that this choice of reference metric f μν , first, cannot produce any infinite value for the bulk action (3.1) (since the bulk action only contains non-negative powers of f μν ), second, it does not preserve general covariance in the transverse spatial coordinates x 1 , x 2 , . . . (since f μν only depends on the spatial components h i j of the spacetime metric). Regarding Eqs. (3.3) and (3.9), the interaction terms U i 's can be calculated as Using the electromagnetic field equation (3.5), the gauge potential is obtained as follows where q is an integration constant related to the electric charge. As a result, the non-zero components of the Faraday tensor are Now, by use of the gravitational field equations (3.4), we are going to obtain the AdS black hole solutions. It is well known that solutions of TOL gravity with different Lovelock coefficients α 2 and α 3 are mathematically too long; therefore not appropriate to be studied. These complete solutions have been proposed in [70]. The black hole solutions can still be found if α 2 and α 3 are dependent to each other. In order to have practical black hole solutions, we consider the following special case for Lovelock coefficients which is a standard simple choice [70][71][72][73][74][75]. Considering Eq. (3.13), one may find one real and two complex (conjugate) solutions for the metric function ψ(r ). The real valued solution for the metric function ψ(r ) is obtained as (3.14) with in which the parameter m 0 is a positive integration constant related to the finite mass of the obtained spacetimes. The above metric function satisfies simultaneously all components of the field equations (3.4) and reduces to the Lovelock-Maxwell black hole solution as m → 0 [70]. In addition, for α → 0, static black hole solutions of (Einstein) massive gravity [40] can be recovered as follows (3.16) The existence of the possible curvature singularity could be explored by use of calculating the Kretschmann scalar which is (3.17) Taking into account the metric function ψ(r ), one may find R αβγ δ R αβγ δ ∝ r −4d 2 near the origin (r → 0). In addition, the Kretschmann scalar is finite for r > 0 and diverges at r = 0, and in conclusion, we regard the origin as an essential singularity of the curvature. This physical singularity can be covered by an event horizon. The roots of the metric function ψ(r ) = g rr = 0 specify the number of horizons. Surprisingly, numerical calculation shows the metric function ψ(r ) could have more than two roots for all horizon topologies in contrast to the usual solutions of Lovelock and Einstein gravities (reported in Refs. [41,44,134] as well). Evidently, this is due to the massive interaction terms. In Fig. 5, we have depicted diverse cases for the possibility of the existence of the multi-horizon solutions in LM gravity. In the left panel of Fig. 5, we have displayed a typical example for the behavior of the metric function ψ(r ) in Lovelock gravity (with zero mass for gravitons). As shown in Fig. 5, in LM gravity (right panel), the metric function may have (i) four (ordinary) roots, (ii) one extreme and two roots, (iii) one extreme root, (iv) two roots, or (v) without any root (which all the roots are real and positive). Hence presented solutions may be inter-preted as the black holes with (three) four horizons, extreme black holes, black holes with one non-extreme horizon, or naked singularity. Hereafter, we assume that r + is the event horizon radius of the black hole solutions (3.14) and can be numerically computed by finding the largest real positive root of ψ(r ) = 0. In order to investigate the asymptotic behavior of spacetime, we consider the metric function ψ(r ) at large r , i.e.
We can find that the obtained solutions are asymptotically We are more interested in studying AdS black holes since these types of black objects admit dual interpretation and also possess certain phase transition(s) in the extended phase space. Hereafter, we assume only AdS black holes.

Thermodynamics of LM charged-AdS black holes
Here, we examine traditional form of the first law of black hole thermodynamics in which the cosmological constant is considered as a fixed parameter in the theory. We calculate conserved charges and thermodynamic quantities associated with the charged-AdS black hole solutions (3.14) in LM grav-ity. First, using the definition of surface gravity, Eq. (2.12), the Hawking temperature can be obtained as where For the black hole entropy, since the area law is generally not valid for higher curvature gravities, one should use another approach to calculate it. For asymptotically flat spacetimes, one can easily apply the Wald method [135] as whereR μνγ δ ,R μν andR are, respectively, the Riemann tensor, the Ricci tensor and the Ricci scalar of the induced metric g μν on (d − 2)-dimensional boundary. The modified entropy of TOL gravity is obtained as  (3.26) where r + = r + (S) and as introduced in Eq. (3.15). Now, it can be straightforwardly checked that by computing the following intensive quantities we obtain the same quantities as those found before for T (3.19) and (3.24). In conclusion, we have a prescription at hand for the first law of black hole thermodynamic which takes the form This prescription for the first law is not consistent with Smarr formula based on scaling argument. In the next section (Sect. 3.4), we reconsider the first law of thermodynamics in the extended phase space and resolve this inconsistency.

Extended phase space thermodynamics
In this section, by treating the negative cosmological constant as a thermodynamic pressure, we reconsider the first law of black hole thermodynamics Eq. (3.30) and makes it consistent with the Smarr formula. In order to extend the first law of thermodynamics for LM gravity, we also regard the massive couplings c i and Lovelock coefficient α as thermodynamic variables.
The thermodynamic quantities M, T , S, Q and have been already derived in Sect. 3.3. The conjugate quantity of the thermodynamic pressure, i.e., the thermodynamic volume, is defined as Now, taking into account c i 's and α as thermodynamic variables, the finite mass M (enthalpy) will be a function of new variables, i.e. M = M(Q, S, P, α, c i ). Consequently, the first law of thermodynamics in the extended phase is rewritten as and C i 's are conjugate quantities corresponding to the thermodynamic variables c i with the following explicit forms Moreover, after a long and tedious calculation, it can be derived that obtained thermodynamic quantities obey the Smarr relation as in which (3.36) Since c 2 is a dimensionless parameter in the massive theory of gravity, it does not appear in the Smarr formula and has no thermodynamic contribution. Furthermore, one can get the Smarr relation by invoking the method of scaling argument and performing the dimensional analysis for all variables in the theory as As a result, we find in which, evidently, c 2 has scaling weight equal to zero and the conjugate potentials k consist of three terms originating respectively from the dependence of the entropy, the mass and the bulk Hamiltonian on the Lovelock couplings α k 's (for more details see [136]). However, the mathematical structure of Regarding the condition (3.13), the general form of Smarr Here, we should warn the readers that, so far, the unusual variations of c i and α originate only from the requirement of the consistency between the first law and the Smarr formula (see Eqs. (3.32) and (3.35)). For example, in the presence of the Born-Infeld electrodynamics, a new thermodynamic quantity B (interpreted as Born-Infeld vacuum polarization) conjugate to the nonlinear electromagnetic parameter β is defined for AdS black holes [93]. But, in the view of AdS/CFT correspondence, more justification is needed to find the correct first law. In this regard, the first law of thermodynamics for general asymptotically locally AdS black hole spacetimes can be proven using Noether's theorem and covariant phase space techniques which is described in detail in Ref. [137]. As first suggested in Ref. [87], the inclusion of as a thermodynamic variable may find application in the AdS/CFT context, and, consequently, variation with respect to on the gravity side essentially amounts to changing the number of colors N on the field theory side (for more detailed discussion see Ref. [138]). In addition, the holographic Smarr relation beyond the large N limit is investigated in Ref. [139] which proves that the bulk correlates of subleading 1/N corrections to the Smarr relation are related to the couplings in Lovelock gravity theories.

Thermal stability of LM charged-AdS black holes
In this section, we perform a thermal stability analysis for the extended phase space of LM charged-AdS black holes in the canonical ensemble, so the conserved (electric) charge Q, Lovelock coefficient α, massive couplings c i and the pressure P will be regarded as fixed thermodynamic quantities. In this ensemble, we will show the heat capacity at constant P, Q, α and c i by the symbol C P,X i , in which X i denotes the extensive quantities Q, α and c i . Since negativity of heat capacity, i.e. C P,X i < 0, indicates local thermodynamic instability, we look for regions where C P,X i is positive. For this purpose, we calculate the heat capacity as where (3.41) As mentioned before, we have restricted our study to U i up to the fourth interaction term (U 4 ). To study local thermodynamic stability of black holes, first, we explore the physical temperature for topological black holes. Numerical calculations for thermal analysis show that for a black hole with definite mass M, there always exists a lower value for the radius of event horizon, i.e. r b , in which the black hole temperature is always positive for r + > r b . According to Eq. (3.19), it can be seen that temperature always has a root regardless of horizon geometry. In Fig. 6, the typical behavior of temperature is depicted for spherical black holes with various values for spacetime dimensions (d), Lovelock coefficient (α) and massive parameter c 0 . As seen, the value of r b depends on many parameters in the theory. For instance, it is an increasing function of spacetime dimensions (d) and a decreasing function of α and c 0 . In addition, we plot Fig. 7 to investigate the qualitative behavior of temperature for Ricci flat and hyperbolic black holes. For Ricci flat black holes, the behavior of temperature is qualitatively similar to spherical black holes (see left panel of Fig. 7). The result is radically different for hyperbolic black holes. In this case, interestingly, an infinity is observed at the divergence point r i = √ α. Again, the black hole temperature is always positive for r + > r b (see middle and right panels of Fig. 7). If r i > r b , an infinite (positive) temperature is observed for hyperbolic black holes with horizon radius of r + = r i = √ α. It should be emphasized that the temperature behaves as T ∝ r + for large values of the event horizon r + (see Eq. 3.19). Therefore, temperature diagrams (Figs. 6 and 7) diverge at r + → ∞ for LM AdS black holes with diverse horizon topologies. Now, we discuss positivity of the heat capacity with respect to the event horizon topology (k) case by case and in detail.
Spherical black holes (k = 1): In Figs. 8, 9, 10, we have displayed various cases for the behavior of C P,X i with respect to r + . As seen, the heat capacity always has only one root (r b ) in which it is negative definite for regions r + < r b . For regions r + > r b , as before in Einstein gravity, there are various possibilities. Depending on the values of q, α, m (graviton mass), massive gravity couplings (c 1 , c 2 , c 3 , c 4 ) and spacetime dimensions (d), the heat capacity (i) may be an increasing function of r + , (ii) may have two divergence points, or iii) has solely one divergence point. LM AdS black holes would be thermally stable if C P,X i were an increasing function of r + without any divergence or root in regions r + > r b . For the second possibility, again, we refer to the first and second divergence points as r m and r u respectively (r b < r m < r u ). As seen in Figs. 8, 9, 10, for regions r b < r + < r m and r + > r u , LM AdS black holes are thermally stable and for regions r m < r + < r u are unstable. The third item (iii) appears in thermodynamic phase space of all kinds of topological black holes, so we postpone this case until later in this section.
For the sake of completeness, Figs. 8, 9, 10 are plotted for different cases to investigate the effects of Lovelock coefficient and massive graviton parameter (m) on thermal stability of the solutions in higher dimensions. Since the quantities m and c 0 are always coupled, thus we would rather vary c 0 instead of m. The results are qualitatively similar. We found that thermally stable regions are drastically affected by both of the deformation parameters, i.e., α and m.
Moreover, since the thermodynamic phase space of topological black holes have been enriched due to the higher order curvatures and the graviton's mass, we can speculate that the heat capacity may have more than two divergence points which indicates more complex critical behaviors with the associated phase transitions. Therefore, by fine tuning   the thermodynamic variables, one can form other possibilities besides the three items which were presented previously. These possibilities can be presented as: item (iv) three diver-gence points for the heat capacity and item (v) existence of four divergence points. In these cases, the absence of any root (i.e. r b in the previous cases) for the heat capac-  Fig. 11 for the corresponding P −r + , T −r + and G − T phase diagrams ity plays an essential role. The case related to the item (iv) is depicted in Fig. 11 which shows that there are two unstable regions and two stable ones denoted by SBHs and LBHs. More investigations, which were also done for another black hole set-ups, reveal that this critical behavior corresponds with the RPT. In Fig. 24 in Sect. 4.4, we have plotted the corresponding P − r + , T − r + and G − T phase diagrams which confirm the mentioned RPT behavior. The behavior associated with the four critical points, item (v), is displayed in Fig. 12. As shown, there are two unstable regions and three stable ones denoted by SBHs, IBHs and LBHs. Thus, the well-known small/intermediate/large black hole phase transition is inferred easily. The corresponding P −r + , T − r + and G − T phase diagrams are plotted in Fig. 25 which support the existence of triple points and the associated small/intermediate/large black hole phase transition as expected. It should be emphasized this behavior, item (v), is only observed in the spherical black holes of LM gravity.  8 and 13), but they are completely different from their counterparts in Einstein gravity, in which the heat capacity cannot possess any divergence point. In fact, depending on the values of parameters in the theory, the heat capacity may possess infinities (at one, two or three points) like those in the case of black holes with spherical horizons (items i, ii, iii and iv). In this regard, according to Fig. 13, black holes are thermally stable in regions r b < r + < r m and r + > r u , and unstable for r m < r + < r u (the same notations as the case of spherical black holes, items i and ii, have been assumed here). In addition, the case of three divergence points can take place the same as the spherical black holes (see Fig. 11) signaling a RPT (the analytic conditions to have the RPT behavior are presented in Sect. 4.3). Hyperbolic black holes (k = −1): Interestingly, in the case of k = −1, there exists two positive roots and one divergence point (r m ) for the heat capacity (not seen in Einstein gravity). It is a general consequence of higher curvature terms of TOL gravity and the results are radically different from hyperbolic black holes in Einstein gravity. As shown in Fig. 14, the sequence of roots and divergence point is highly important to find stable regions. The smaller and larger roots are referred to as r b 1 and r b 2 , respectively. Investigations show that two sequences are possible: (i) the divergence point is smaller than the roots (r m < r b 1 < r b 2 ), and therefore, LM AdS black holes are thermally stable only in regions r m < r + < r b 1 and r + > r b 2 ; (ii) the divergence point is larger than the roots (r b 1 < r b 2 < r m ), and hence, black hole solutions are thermally stable only in regions r b 1 < r + < r b 2 and r + > r m . In the other regions, hyperbolic black holes behave thermally unstable. The existence of two divergence points for the heat capacity of hyperbolic black holes will be clarified in a moment. Now, we seek other qualitative behaviors of the topological black holes' heat capacity. As mentioned in Sect. 2.3, the heat capacity of spherically symmetric AdS black holes in Einstein gravity may have solely one divergence point (r m ) which is positive around such divergency. Interestingly, that behavior can occur for the obtained black hole solutions with various horizon topologies in LM gravity. In Fig. 15, this possibility is depicted for all topological black holes in 7-dimensions. Such single divergence point, with positive C P,X i around it, may signal critical behavior of the topological black holes. As seen, all plots have the same qualitative behavior as the spherically symmetric AdS ones in Einstein gravity. The only difference is the asymptotic behavior of the heat capacity of the hyperbolic black hole in comparison with spherical and Ricci flat black holes. In fact, when r + → ∞, C P,X i approaches large negative values for hyperbolic black hole solutions. Instead, large positive values are observed for the asymptotic behavior of the heat capacity of AdS black holes with spherical and Ricci flat horizons. Consequently, there is a lower bound, referred to as r b (r b < r m ), where the heat capacity is positive for spherical and Ricci flat black holes in regions r + > r b . For hyperbolic black holes, there is an upper bound for the horizon radius (r u ) besides the lower bound (r b ), in which the heat capacity is positive in regions r b < r + < r u and negative for the other regions. In addition, there is another possibility for the case of one divergence point in the heat capacity, which indicates the Hawking-Page phase transition. According to Fig. 16, around such divergency, the heat capacity is negative and positive corresponding to the area on the right and the area on the left, respectively. In all cases, besides a divergence point (r m ), C P,X i may have one or two positive roots (referred to as r b 1 and r b 2 ) depending on the chosen parameters. When there exist two roots, the divergence point is located between the roots. This situation is different from the previous case of hyperbolic black holes (see the items i and ii related to hyperbolic black holes). Assuming that there are two roots (r b 1 < r b 2 ), an unstable region is observed for r m < r b 2 in the heat capacity's plots of Fig. 16. Consequently, thermally stable regions correspond with r b 1 < r + < r m and r + > r b 2 .

Critical behavior and vdW phase transition
In this section, we study the critical behavior of the LM charged-AdS black holes in the canonical ensemble under certain conditions which vdW phase transition appears. Geo- and r + is a function of thermodynamic volume, V . This is a worthwhile equation since one can easily recover the equations of state of charged-AdS black holes in Einstein, Gauss-Bonnet and massive gravities (and also any possible combination of them). It is well-known that black hole solutions in TOL gravity have no Gauss-Bonnet limit [140], thereby one cannot obtain Gauss-Bonnet gravity's outcomes by taking limit α 3 ≡ α TOL → 0 in the context of TOL gravity. But there is a tricky procedure, for example, to obtain the equation of state of AdS black holes in Gauss-Bonnet gravity from TOL gravity. Since the terms including α 2 are due to the effect of third order term in the Lovelock Lagrangian (3.2), one can recover the Gauss-Bonnet-massive equation of state for charged-AdS black holes in Ref. [45] by taking limits α → d 3 d 4 α GB and α 2 → 0 according to the Lovelock coefficient conditions (3.13) which state that α 2 ≡ α GB ∝ α and α 3 ≡ α TOL ∝ α 2 . It should be noted that in the above discussion we did not assume the square of Gauss-Bonnet coefficient (α GB 2 ) is too small. In what follows, we restrict our study to U i up to the fourth interaction term (U 4 ). Again, the physical equation of state can be obtained by translating the geometric version, Eq. (4.1). The same result is obtained for associated specific volume (by performing the steps as stated in 2.4). To do that, one has to define the shifted Hawking temperature [107][108][109] as Rewriting the geometric equation of state (4.1) as k BTphys + · · · . (4.6) Hence, we can identify the specific volume v as v = 4r + d 2 P /d 2 . Since in the geometric units r + = d 2 v/4, hereafter, we would rather use the event horizon radius (r + ) instead of the specific volume (v) in order to analyze the criticality. However, there are physical differences between r + , v and V , but because of their dependency to each other, criticality in one of P − r + , P − v or P − V planes indicates criticality in the others.
Let's compare the LM equation of state (4.4) with the equation of state of the RN-AdS black holes (2.29). In the right hand side of the RN-AdS equation of state, the first and the last terms are always positive and the second term can be positive, zero or negative. In fact, signs of the presented terms in the equation of state could ensure the critical behavior for a given black hole solution. We introduce the signature of the equation of state of RN-AdS black holes as P(+, ±, +). We saw that there does not exist P − V (P − r + ) criticality for black holes with Ricci flat and hyperbolic horizons. In order to have P − V (P − r + ) criticality two positive terms and one negative term are needed, i.e., an equation of state with signature P(+, −, +) which is the case with spherical horizon (k = +1). Therefore, at least, two positive terms and one negative term in the equation of state can possibly ensure the critical behavior and phase transition for a given black hole. Regarding this, the equation of state of LM charged-AdS black holes (4.4) with signature P(±, ±, ±, ±, +, ±, +) predicts critical behavior and phase transition for Ricci flat and hyperbolic black holes as well as spherically symmetric black holes depending on the massive coupling coefficients (c i ), Lovelock coefficient (α) and topological factor (k).
We are looking for the inflection point of isothermal P−r + diagrams, the subcritical isobar of T − r + plots, and the characteristic swallow-tail form of G − T diagrams for the obtained black hole solutions according to Sect. 2.4. These pieces of evidence guarantee the existence of phase transition and indicate the vdW like behavior for the LM AdS black holes. In our considerations, we suppose that all the massive coupling coefficients are simultaneously positive (c 1 = c 2 = c 3 = c 4 = +1) or negative (c 1 = c 2 = c 3 = c 4 = −1) and keep track of the effect of higher order curvature terms of the Lovelock Lagrangian on the outcomes of massive gravity. Later, we do not impose this assumption and will summarize the results of arbitrary signs for the massive couplings (c i = ±1).
First, we study the P−r + diagrams of LM AdS black holes with various topological factors. In an isothermal P − r + diagram, the critical point is an inflection point and can be obtained from Eq. (2.35). Investigation of the critical behavior is not possible, analytically, and so we apply the numerical analysis. Using Eqs. (4.1) and (2.35), the critical point can be found numerically. At the critical point, we refer the value of horizon radius as a critical horizon radius, r c . The critical values for the topological black holes in d = 7 dimensions are gathered in Table 1 in which we have focused on the effect of horizon topology factor (k) while keeping other parameters fixed.
Considering the obtained critical values in Table 1 for the pressure, horizon radius and temperature, one can plot their corresponding phase diagrams. In the left panels of Figs. 17,18,19, the characteristic behavior of pressure as a function of event horizon radius (r + ) is depicted for the topological black holes. Compared to P −r + diagram of vdW fluid or RN-AdS black holes (see Fig. 4), it is seen that the associated P − r + diagrams for LM AdS black holes qualitatively behave like vdW fluid. Therefore, critical radius (inflection point) can be found for Ricci flat or hyperbolic black holes as well as spherically symmetric black holes. For all P − r + diagrams, the temperature of isotherms decreases from top to bottom. For T > T c , the isotherms correspond to the ideal gas with a single phase. For T < T c , a two-phases behavior is seen, and in comparison with the liquid/gas system, there exists (first order) small/large black hole phase transition for topological black holes. It is notable that, based on Maxwell's equalarea law, the (unphysical) oscillatory part of each isotherm is replaced by a line of constant pressure. Now, we are looking for the characteristic swallow-tail form of G − T diagrams. The Gibbs free energy is found by computing on-shell action (I G ) or using the Legendre transformation, G = M − T S. Analytical calculation of the Gibbs free energy is too lengthy to write here and, therefore, we leave out the analytical result for reasons of economy. We have plotted the Gibbs free energy as a function of temperature for various pressures in right panels of Figs. 17, 18, 19. As seen, obviously, G − T diagrams indicate the characteristic swallow-tail behavior for all types of topological For the sake of completeness, the qualitative behavior of temperature as a function of horizon radius (which corresponds to specific volume, v) are depicted in the middle panels of Figs. 17, 18, 19. Comparing Fig. 4 with Figs. 17, 18, 19, T − r + diagrams shows a vdW like behavior. Evidently, Ricci flat and hyperbolic black holes can qualitatively imitate the critical behavior of spherical black holes in LM gravity. It should be emphasized that the vdW like behavior persists in higher dimensions for all kinds of topological black holes. In addition, the universal ratio (i.e., P c r c T c ) is a function of event horizon topology.
To be more specific, we analyze the equations of state and phase transitions for topological black holes case by case and in details. We focus on the effects of Lovelock coefficient (α), graviton mass parameter m, spacetime dimension (d) and topological factor (k) and present related tables (see Tables 2, 3 , 4, 5, 6, 7, 8, 9 and 10). In this regard, we reveal a peculiar phase transition and critical behavior for hyperbolic black holes in the LM gravity in higher dimensions of spacetime (i.e., d 7).

Spherical horizon(k = +1):
The LM equation of state for spherical black holes with signature P(±, ±, ±, ±, +, −, +) reads (4.7) Our investigations show this equation of state predicts the critical behavior in higher dimensions. In this case, if all the massive coupling coefficients (c i ) be positive (negative) definite, one (two) physical critical point(s) can be found at most. In Tables 2, 3 and 4, the critical values for pressure, horizon radius and temperature have been computed for various spacetime dimensions (d), Lovelock coefficient (α) and graviton mass parameter (m). According to Tables 2 and 4, critical pressure and temperature are increasing functions of d and m whereas the critical horizon radius is a decreasing function of them. The universal ratio P c r c T c is an increasing function of spacetime dimensions the same as RN-AdS black holes in Einstein gravity (see Eq. 2.37).
In addition, the functional form of critical values with respect to the Lovelock coefficient (α) are investigated in Table 3. According to this table, critical pressure and temperature are decreasing functions of α and critical horizon radius is an increasing function of it. Now, suppose there is a certain vdW phase transition for a given spherical black hole   Table 3 Spherical black holes: k = 1, d = 7, q = 1, m = 1, and  set-up in the absence of higher curvature terms, α = 0 (this case has been illustrated in the first line in Table 3). Interestingly, in the presence of higher curvatures (α = 0), criticality is observed for all values of α. We examined the asymptotic behavior of critical values with respect to α (i.e., α → ∞) and found that critical pressure, temperature and universal ratio approach zero, and critical radius moves towards infinity (i.e., r c → ∞). These results have been confirmed in higher dimensions as well. We will see that this is not the case for hyperbolic black holes in LM gravity.
Ricci flat horizon (k = 0): The LM equation of state for Ricci flat black holes is given as  with the signature P(±, ±, ±, ±, +). As seen in Eq. (4.8), the effect of higher order curvatures of TOL gravity, which encodes in the Lovelock coefficient α, vanishes for Ricci flat black holes, and the only effect of them comes from the location of the event horizon according to ψ(r + ) = 0 (see Eq. 3.14). In this case, one can show that there exists only one critical point depending on all massive coupling coefficients be positive or all of them be negative definite (see Sect. 4.3). Tables 5 and 6 show some Ricci flat black hole set-ups with various values of d and m. Interestingly, as shown in Table 6, there is a lower value for the graviton mass parameter, referred to as m b , in which no phase transition takes place in region m < m b . Remarkably, Ricci flat black holes can experience critical behavior and small/large black hole phase transition by giving mass to gravitons.

Hyperbolic horizon (k = −1):
The LM equation of state for hyperbolic black holes reads with the signature P(±, ±, ±, ±, +, +, +).Amazingly, when all the massive coupling coefficients are positive or all of them are negative, two (physical) critical points can be found at most. In fact, the equation of state (4.9) could have three positive roots, but only two of them can be physical. We refer to the smaller and larger critical horizon radii as r c 1 and r c 2 respectively (r c 1 < r c 2 ). It should be noted that the third critical point, which actually has the smallest value for the critical horizon radius, is always unphysical since the associated black hole has a negative temperature. In the Tables 7 and 8, the smaller and larger critical horizon radii (r c 1 and r c 2 ) have been computed for the LM AdS black holes with hyperbolic horizons in higher dimensions. Hyperbolic black holes with smaller critical horizons (r c 1 ) experience large values for the temperature and pressure, and therefore, we call them as "hot black holes". In this sense, hyperbolic black holes corresponding to larger critical horizons (with small values for the temperature and pressure respect to the hot black holes) are referred to as "cold black holes". In Figs. 20 and 21, we have depicted the typical behavior of P − r + , T − r + and G − T curves in the vicinity of the first and second (physical) critical points corresponding to r c 1 and r c 2 respectively. Also, the associated critical  Tables 7 and 8. In Fig. 20, for the first critical point (corresponding to the smaller horizon radius, denoted by P c 1 , r c 1 and T c 1 ), we observe the characteristic swallowtail behavior in G − T diagrams for pressures in the range P > P c 1 in contrast to the vdW phase transition which only takes place for P < P c . From this diagram, it can be inferred that a first order phase transition occurs for P > P c 1 since the first derivative of the Gibbs free energy at the critical point is undefined. In P −r + diagrams, interestingly, the (unphysical) oscillating part of isotherms takes place for T > T c 1 which means the existence of two-phases behavior, and the oscillating part is replaced by an isobar according to Maxwell's equal area law. For region T < T c 1 in P − r + diagrams, the one phase behavior corresponding to ideal gas is observed. Moreover, T − r + plots reveal that oscillating part of critical isobars occurs for the pressures in the range P > P c 1 .
Comparing with the vdW phase transition, this critical behavior is completely reverse. This evidence shows hyperbolic black holes could potentially experience the reverse vdW   like behavior for (first order) phase transition at high temperature and pressure which is a remarkable result. Further, in the next part of this section, we will uncover the theory dependency's origin of the "reverse behavior" in LM gravity with a detailed discussion. As far as we know, there is no reverse vdW phase transition in usual thermodynamic systems. In Fig. 21, the qualitative behavior of the hyperbolic (cold) black hole at the second critical point (r c 2 ) is displayed. At this critical point, we observe the standard vdW phase transition which explained in details before. Interestingly, numerical calculations, which are presented in Tables 7 and 8, show that the critical pressure, horizon radius, temperature and the universal ratio ( In addition, according to numerical calculations which are presented in Table 10, there is an upper limit for the value of Lovelock parameter, α u , in which no phase transition could happen for α > α u . This statement does not hold for LM AdS black holes with spherical horizon which is inferred from Table 3. In Table 10, we have shown only those branches corresponding to standard vdW phase transitions and the associated larger critical radii. As α → α u , the two critical points associated with vdW and reverse-vdW phase transitions move towards each other, and eventually, when α = α u both phase transitions with the associated critical radii disappear. This shows that inclusion of higher curvature terms (based on Lovelock Lagrangian) affects drastically the criticality of AdS black holes with the hyperbolic horizon in massive gravity; by tuning the Lovelock coefficient (α) the first-order phase transition can be produced or ruined (see Table 10). Moreover, the universal ratio P c r c T c is a decreasing function of α in the range 0 < α < α u .
Until now we have considered the equations of state of AdS black holes in the LM gravity framework and disclosed a peculiar critical behavior and a strange phase transition for hyperbolic black holes at high temperatures. In order to have a better understanding of the theory dependency and nature of phase transitions, we will separately study the equations of state of AdS black holes in Lovelock and massive gravities.
Here, we confront with an interesting situation: for spacetimes with d > 12, there does not exist any critical point for Lovelock AdS black holes with spherical horizon geometry, and so phase transition does not take place. In 7-dimensions, we observe only one critical point for spherical black holes, and in 8 d 11 there are always two critical horizon radii, r c 1 and r c 2 , which the corresponding critical points can be physical or unphysical (the critical data are presented in [97]). As a result, in 7-dimensions, the vdW behavior, and in 8 d 11, the reentrant behavior for phase transition are observed. In d = 12 we do not observe P − V (P − r + ) criticality and instead, there is a cusp-like behavior in the G − T diagram. The effect of the U (1) charge could drastically alter the number of critical point(s) and the nature of phase transitions. In fact, in spacetime dimensions with the range 8 d 11, there is a lower value for the electric charge (Q b ), where for Q > Q b , one of the critical points disappears and consequently the reentrant behavior is replaced by the vdW like phase transition (some critical data associated with the Lovelock charged-AdS black holes in higher dimensions are presented in [99]). Interestingly, for d 12, one critical point emerges when the U (1) charge is considered and the vdW phase transition takes place. In addition, in d = 7, the inclusion of the U (1) charge only changes the location of the critical point. In Fig. 22, as an example, we have plotted the critical behavior of a spherical black hole. As seen, in this case, criticality in P − r + plane is qualitatively similar to the critical behavior of the vdW fluid, RN-AdS and LM AdS black holes.
In the case of hyperbolic black holes, there always exists only one critical radius (r c = √ α) in all spacetime dimensions (d 7) according to Eq. (4.12). As stated in Sect. 3.5, the temperature of hyperbolic black holes blows up at the point r + = r i = √ α which is referred to as the thermodynamic singularity [100] since all isotherms cross at r c = √ α = d 2 v c /4. Also, the heat capacity is zero at this point. The critical point corresponding to this thermodynamic singularity is called the isolated critical point and is regarded as the first example of a critical point with non-standard critical exponents as α = 0, β = 1, γ = 2 and δ = 3 [100,101] which are different from those of vdW fluid with α = 0, β = 1/2, γ = 1 and δ = 3. Since critical exponents determine the behavior of thermodynamic quantities near the critical point, so various critical exponents imply different behaviors in phase diagrams. As a result, in the case of hyperbolic black holes, a stretched swallow-tail behavior (also referred to as cusp-like behavior) in G − T diagram is observed for pressures in the range P > P c , in which P c is the pressure at the thermodynamic singularity. When the U (1) charge is considered, one additional critical point might emerge. Regarding the charged case, if the value of U (1) charge is more than a lower value (Q > Q b ), there exist two critical radii which the smaller one is always unphysical (the corresponding black hole has a negative value for the temperature) and the larger one can be physical. For the small enough charges (i.e., 0 < Q < Q b ), those critical points are created near the isolated critical point, and the both of them could be physical (r c 1 and r c 2 ). We observe that by increasing the value of the electric charge (Q) the distance between the thermodynamic coordinates of those points is increased in the extended phase space. As a matter of fact, one can grow up the thermodynamic quantities associated to the critical point by increasing the U (1) charge and produce a first order phase transition at high temperature and pressure for hyperbolic black holes.
Here, an interesting phenomenon emerges and persists in higher dimensions (d ≥ 7). Our investigations show that the reverse vdW phase transition is a characteristic feature of Lovelock AdS black holes with hyperbolic horizon corresponding to those larger critical points (this strange behavior is already pointed out in [97,98,100]). As seen, Fig. 23 exposes the origin of the reverse vdW like behavior which has been found in Sect. 4.1 for hyperbolic black holes in the LM gravity. In Fig. 23, we observe the existence of inflection point in the isothermal P − r + diagrams, the subcritical isobar of T − r + plots, and the characteristic swallow-tail form of G − T diagrams, but, in contrast to the behavior of vdW fluid, in the opposite way. Therefore, we conclude higher order curvature terms based on the Lovelock Lagrangian are responsible for the reverse vdW phase transition. It should be emphasized that the Lovelock equation of state (4.10) has been obtained by the assumption of the Lovelock coefficient condition (3.13). In the more general case where Lovelock coefficients are independent, one may obtain three critical points for black holes which indicate the appearance of triple point [100].

Massive gravity: the phase transition revisited
In the context of massive gravity, the equation of state of charged-AdS black holes is given as in which we have used the zero higher curvature limit (α → 0) of the LM equation of state (4.4) and the shifted Hawking temperature is the same as before, i.e., T = T − m 2 c 0 c 1 /4π . As seen, the topological factor (k) in the RN-AdS equation of state (2.29) is replaced by the combination (k + m 2 c 2 0 c 2 ) which suggests possible critical behavior for black holes with non-trivial horizon topologies. Clearly, the massive charged-AdS equation of state (4.13) with signature P(±, ±, ±, ±, +) admits phase transition for black holes with spherical, Ricci flat and hyperbolic horizon geometries [106]. In 4-dimensional spacetime, the critical point is obtained from the positive root of the equation r 2 + (k + m 2 c 2 0 c 2 ) − 6q 2 = 0, in which we applied Eqs. (2.35) and (4.13). Hence, there exists a critical horizon radius (r c ) when (k + m 2 c 2 0 c 2 ) 0. It should be noted that, like RN-AdS black holes in Einstein gravity, the U (1) charge is necessary to have critical behavior and phase transition in a 4dimensional spacetime. On the other hand, as indicated in [107], in higher dimensional spacetimes (d 5), the U (1) charge is unnecessary for the appearance of criticality in massive AdS black holes since the third and the fourth massive potential terms in the right hand side of Eq. Interestingly, the critical point(s) is independent of the first massive coupling coefficient (c 1 ). Considering the uncharged case (q = 0) for the sake of simplicity, the critical point equation (4.14) can be solved simply as .  .18), when all the massive coupling coefficients are simultaneously positive (negative) definite, there exists one critical radius and can be determined using Eq. (4.17). In order to have two critical points, one should consider some specific signs for the massive couplings (c i ) based on Eq. (4.18). When the combination (k + m 2 c 2 0 c 2 ) is positive, two critical points can be found assuming that c 3 < 0 and c 4 > 0, and when (k + m 2 c 2 0 c 2 ) < 0, one has to assume c 3 > 0 and c 4 < 0. In the both cases, the massive coupling coefficients must satisfy the constraint 3d 4 m 2 c 2 0 c 2 3 > 8d 5 c 4 (k + m 2 c 2 0 c 2 ). The existence of two critical points for (neutral) AdS black holes is possible when the spacetime dimensions are more than five (d 6) and an interesting phenomenon emerges which is called the RPT (for more details see [108]). The inclusion of nonlinear electromagnetic fields can increase the number of critical point(s) [93], and as a result, in the context of Born-Infeld-massive gravity [109], the RPT appears in 4-dimensions and the so-called triple point in spacetime dimensions more than five (d 5). It should be noted these considerations were done in the canonical ensemble.

LM gravity: RPTs and triple points
This section is devoted to studying the possibility of the appearance of the RPT and triple point in the phase structure of the LM AdS black holes. In the previous section, we indicated that under certain conditions the equations of state of topological black holes in pure massive gravity (without non- and T > T c1 (dashed line). Middle panel: P c2 < P < P c1 (continuous lines), P = P c1 (dotted lines) and P > P c1 (dashed lines). Right panel: P tr < P < P z (continuous lines), P = P c1 (dotted lines) and P > P c1 (dashed lines). Critical data: (P c2 = 2.74376, r c2 = 0.131663, T c2 = 1.74384), (P c1 = 4.47518, r c1 = 0.249414, T c1 = 1.83107), (P tr 3.93852, T tr 1.76669), and (P z 3.96314, T z 1.76709). Note: The values of Gibbs free energies in the G − T diagram (continuous and dashed lines) have been slightly shifted up to avoid isobaric curves overlap with each other. The corresponding behavior of the heat capacity near the virtual triple point (P tr ) is depicted in Fig. 11 trivial electromagnetic fields like BI electrodynamics) may have up to two critical points and thus exhibit vdW and RPTs which the latter corresponds to three-phases behavior. In Sect. 4.2, we showed Lovelock (un)charged-AdS black holes with spherical horizon can exhibit vdW and RPTs. Also, for the Lovelock (un)charged-AdS black holes with the hyperbolic horizon, the reverse vdW like behavior is observed which can be accompanied by a (normal) vdW phase transition. Consequently, since there are many thermodynamic variables in the extended phase space of the LM AdS black holes, one expects these black holes may enjoy a vast range of thermodynamic behaviors which found in the other gravitational theories similar to those in usual thermodynamics.
Here we intend to examine these possibilities.
First, we consider Ricci flat black holes (k = 0) in the LM gravity. In this case, the effect of higher curvature terms (encoded in the Lovelock coefficient, α) vanishes since α is always coupled to the topological factor k. As a result, using Eq. (4.18), one can find the RPT in spacetime dimensions d 6 for the neutral black holes [108].
Investigation shows the LM AdS black holes with spherical horizon may have up to three physical critical points for charged and uncharged cases. In order to observe the reentrant behavior of phase transition for spherical black holes, we have adjusted the massive coupling coefficients to produce two critical points (referred as r c 1 and r c 2 ) according to Eq. (2.35) and plotted P − r + , T − r + and G − T diagrams in Fig. 24. As a result, a virtual triple point (T tr , P tr ) and another critical point (T z ,P z ) emerge in the phase space. In addition, the behavior of the heat capacity associated to the reentrant and small/intermediate/large phase transitions are depicted in Figs. 11 and 12 near the (virtual) triple point. As seen, by monotonic decreasing the temperature, the black hole system undergoes a RPT for the certain range of pressure (P tr < P < P z ). By another tuning, we arrive at one triple point (T tr , P tr ) and two physical critical points (associated with r c 1 and r c 2 ). This situation is depicted in Fig. 25 in which the Gibbs free energy is displayed near the critical points for various pressures. It should be mentioned that there is a lower value for the U (1) charge (Q b ), where for Q > Q b , one of the critical points disappears. Hence, in the case of charged-AdS black holes, the analogue of the triple point and solid/liquid/gas phase transition can be found only for small enough values of the electric charge, Q.
The phase structure of hyperbolic black holes is really rich and drastically different from those of spherical and Ricci flat horizons. In both charged and uncharged cases, three (physical) critical points can be found for hyperbolic black holes. Interestingly, the analogue of triple point does not exist in the phase structure of these black holes. In fact, besides the existence of the two critical points corresponding with two distinct first order transition, we arrive at an additional reverse vdW phase transition associated to the third critical point in the phase space. This situation is illustrated in Fig. 26 which is a generic feature of this model and persists in all dimensions. This is the first example of such phase structure which is not possible for spherical and Ricci flat black holes.

Concluding remarks
The effects of massive and Lovelock gravities are encoded in the deformation parameters m and α, respectively. In Lovelock massive (LM) gravity, one can simply recover the outcomes of Einstein (by α, m −→ 0), Lovelock (by m −→ 0) and massive (by α −→ 0) theories of gravity. Considering LM gravity, in this paper, we introduced topological black hole solutions and then analyzed thermodynamic properties and critical behavior of AdS black holes in the extended phase space. The asymptotic behavior of the black hole solutions may be (A)dS or flat, and by computing the thermodynamic quantities, we have shown they satisfy the first law of thermodynamics.
Next, by treating the cosmological constant as a thermodynamic variable (pressure), we extended the thermodynamic phase space and proved the massive coupling and Lovelock coefficients, as well as cosmological constant, are required for consistency of the extended first law of thermodynamics with the Smarr formula. In addition, we examined thermal stability in the extended phase space thermodynamics of the LM AdS black holes in the canonical ensemble (where the quantities P, Q, c i and α are held fixed), and showed the qualitative behavior of heat capacity for AdS black holes with different horizon topologies. In this regard, we mainly focused on the topology of event horizons and showed in what regions the topological black holes are thermally stable.
In LM gravity, critical behavior and phase transition occur for all types of AdS black holes (with spherical, Ricci flat and hyperbolic topologies for event horizon) in contrast to Einstein gravity which only admits phase transition for spherically symmetric ones. For Ricci flat black holes, phase transition originated only from the interacting terms of massive gravitons and the effect of higher curvature terms vanishes. Interestingly, we found that there is a lower value for the graviton mass parameter, referred to as m b , in which no phase transition takes place in the region m < m b . This is one of the remarkable characteristics of massive gravity. For hyperbolic black holes, two radically different first order transitions are observed: (i) a (normal) vdW like behavior, and (ii) reverse vdW like behavior. The reverse behavior of vdW phase transition completely comes from the higher curvature terms of Lovelock Lagrangian which is not seen in Gauss-Bonnet gravity (as indicated in [100,141], Gauss-Bonnet black holes with hyperbolic horizon do not admit physical phase transition). The reverse behavior predicts that hyperbolic black holes could experience first order phase transition at high temperature and pressure, which is a novel effect. Moreover, it was shown that the inclusion of higher curvature terms (based on Lovelock Lagrangian) affects the criticality. In fact, for LM AdS black holes with hyperbolic horizon topology, depending on the chosen parameters, there is always an upper limit for the value of Lovelock coefficient (α u ) in which no phase transition could happen for α > α u . In the case of spherical black holes, this statement no longer holds and we observe criticality in the range 0 < α < ∞. Considering tables, we found that the universal ratio, i.e. P c v c T c , is a function of spacetime dimensions (d), topological factor (k), graviton mass parameter (m) and strength of higher curvature terms (captured with Lovelock coefficient, α).
In addition, the vdW, reentrant and analogue of solid/liquid/ gas phase transitions were found in the extended phase space of (un)charged-AdS black holes with the spherical horizon. But, in the case of hyperbolic black holes, reentrant and small/intermediate/large phase transitions were not found. Indeed, the reverse vdW phase transition in the phase space of hyperbolic black holes is accompanied with one or two distinct (standard) vdW phase transitions. To our knowledge, this is the first example of such a phase structure. These pieces of evidence show that the generic features of different theories of gravitation can be summed into a unique model to produce more complex structures for the thermodynamic phase space of black holes.