New perspective for black hole thermodynamics in Gauss-Bonnet-Born-Infeld massive gravity

Following earlier study regarding Einstein-Gauss-Bonnet-massive black holes in the presence of Born-Infeld nonlinear electromagnetic field [S. H. Hendi, B. Eslam Panah and S. Panahiyan, arXiv:1510.00108], we study thermodynamical structure and critical behavior of these black holes through various methods in this paper. Geometrical thermodynamics is employed to give a picture regarding phase transition of these black holes. Next, a new method is used to derive critical pressure and horizon radius of these black holes. In addition, Maxwell equal area law is employed to study the Van der Waals like behavior of these black holes. Moreover, the critical exponents are calculated and by using Ehrenfest equations, the type of the phase transitions are determined.


I. INTRODUCTION
Black hole solution is one of the interesting consequences of general relativity. Although the existence of black holes is vivid, it is an open question to realize interior nature of them in quantitative detail; the main reason comes from the fact that a perfect theory of quantum gravity does not yet exist. Studying the semiclassical phase structure of black holes provides at least preliminary steps for understanding the quantum gravity.
The phase transition plays an important role for exploring the critical behavior of the system near critical point. After the discovery of a phase transition by Hawking and Page [2], black hole phase transitions have been of great interest. It is known that the asymptotically flat vacuum black hole solutions are thermally unstable [3]. While asymptotically AdS black holes are the famous examples of the Hawking-Page phase transition [2] between two stable phases. In order to characterize the critical behavior of a system during the phase transition, one may calculate its critical exponents which are not completely independent. Meanwhile two systems belong to the same universality class if their critical behavior is expressed by the same critical exponents.
The semiclassical phase transition which occurs in the asymptotically AdS spacetimes can be translated to confinement/deconfinement phase transition in the context of AdS/CFT [4]. Regarding the applications of AdS/CFT correspondence in recent years, the similarities between the phase transition of black holes and holographic superconductivity, have achieved a great deal of attention [5,6].
The local (thermal) stability of a system is concerned with how the system responds to small fluctuations of its thermodynamic coordinates. There are various methods that one can employ to investigate the phase structure of a black hole system near its critical point. One of the well-known standard analysis of the locally stability is based on canonical ensemble by studying the specific heats. Phase structure may be also explained by critical quantities that are extracted in the extended phase space. In addition, one may apply the geometrical thermodynamics method for studying the phase transition.
One of the methods for constructing phase structure of a thermodynamical system is through the use of geometry. Meaning, by considering a thermodynamical potential and its corresponding extensive parameters, it is possible to introduce a metric which describes thermodynamical properties of the system. The information regarding thermodynamical properties of the system is extracted from Ricci scalar of the metric. The divergencies of Ricci scalar of the thermodynamical metric mark two important points of the thermodynamical system; bound point and phase transition. There are several methods regarding thermodynamical geometry which are; Weinhold [7,8], Ruppeiner [9,10], Quevedo [11,12] and HPEM [13][14][15]. The geometrical thermodynamics has been employed in the context of different types of black holes [16][17][18][19][20][21][22]. In addition, this method was also used to study phase transition of superconductors [23]. A comparative study regarding different geometrical thermodynamical metrics is done in Ref. [24]. A successful method of the geometrical thermodynamics include all bound and phase transition points in its Ricci scalar through the divergencies. In other words, divergencies of the Ricci scalar and mentioned points must coincide with each other. A mismatch and extra divergency indicate the existence of anomaly which contradict with principles of thermodynamics. Such anomaly was reported for Weinhold, Ruppeiner and Quevedo metrics for different types of black holes [13][14][15]25]. To overcome such problem, HPEM metric was introduced [13]. In this paper, we will regard HPEM method for studying geometrical thermodynamics of black holes under consideration.
Einstein gravity introduces gravitons as massless particles, whereas there are several arguments that state gravitons should be massive particles. In order to have massive graviton, theory of general relativity should be modified to include mass terms. The first attempt for constructing a massive theory was referred to the works of Fierz and Pauli [26] which was done in the context of linear theory. This theory has specific problem which is known as van Dam, Veltman and Zakharov discontinuity. Meaning that propagator of the massive gravity in limit m = 0 is not consistent with one derived for massless case. The resolution to this problem was Vainshtein mechanism which requires the system to be considered in nonlinear regime. In other words, according to the Vainshtein mechanism [27], at some distance below the so-called Vainshtein radius, the linear regime breaks down and the theory enters into a nonlinear framework. Based on this mechanism, the usual general relativity can be recovered from high curvature spacetimes which are introduced with a wide class of non-Einsteinian theories. In the context of a static and spherically symmetric base space, it is also shown that the Vainshtein mechanism can work, correctly, both inside and outside the compact objects [28,29] (See Refs. [30][31][32][33][34] for more details regarding Vainshtein mechanism). On the other hand, generalization of the Fierz and Pauli massive theory to nonlinear one leads to existence of a Boulware-Deser ghost [35]. While solutions to these problems had existed for some time in three dimensional spacetime [36,37], they were not solved in four and higher dimensions. In order to solve such problems de Rham, Gabadadze and Tolley (dRGT) proposed another class of massive gravity [38,39]. Contrary to previous theories, dRGT theory is valid in higher dimensions and it was shown that such theory enjoys absence of the Boulware-Deser ghost [40,41]. This theory build up massive terms by employing a reference metric. A modification in reference metric could lead to another dRGT like massive theory [42]. Black hole solutions of dRGT massive gravity and their thermodynamical properties have been investigated for d-dimensions (d ≥ 3) in Refs. [25,[43][44][45][46][47][48].
On the other hand, one of the well-known theories of higher derivative gravity is Lovelock theory which is a natural generalization of Einstein gravity in higher dimensions. Taking into account the first additional term of Einstein gravity in the context of Lovelock theory (Gauss-Bonnet (GB) gravity), it is believed that GB gravity can solve some of the shortcomings of Einstein gravity [49][50][51]. In addition, GB gravity consists curvature-squared terms which, interestingly, is free of ghosts and the corresponding field equations contain no more than second derivatives of the metric (see Refs. [52][53][54][55][56][57] for more details). Another interesting aspect of GB gravity is that it can be arisen from the low-energy limit of heterotic string theory [58][59][60][61]. Considering GB gravity context, black hole solutions and their interesting behavior have been investigated in many literatures [62][63][64][65][66][67][68][69][70].
On the other hand, one of the main problems of Maxwell's electromagnetic field theory for a point-like charge is that there is a singularity at the charge position and therefore, it has infinite self energy. In order to remove this self energy, in classical electrodynamics, Born and Infeld introduced a nonlinear electromagnetic field [71], with main motivation of solving the infinite self energy problem by imposing a maximum strength of the electromagnetic field. Motivated by the interesting results mentioned above, we study thermodynamic behavior of black holes in GB-massive gravity in the present Born-Infeld (BI) source.
In order to have a better description regarding physics governing a system, it is necessary to decrease different shortcomings of different theories as much as it is possible. This indicates that we should apply more generalizations to solve different shortcomings of theories describing the nature of system. Here, we have considered three generalizations; BI generalization to remove shortcomings of the Maxwell theory, GB gravity to solve different problems of Einstein theory such as renormalization problem, and massive gravity to solve the massless gravitons in both Einstein gravity and GB theory. Such considerations solve some of the shortcomings of theories under consideration, but they also modify the physical properties of the system. In this paper, we intend to investigate these modifications in the context of critical behavior of black holes which in turn provides a reasonable framework for conducting studies in other aspects of physics such as gauge/gravity duality.
The outline of the paper will be as follow. In next section, we introduce action and basic equations related to GB-BI-massive gravity. We also present a brief discussion regarding to the black hole solutions and conserved and thermodynamics quantities. Section III is devoted to study the phase transition through geometrical thermodynamics. In Sec. IV, we investigate critical behavior of the system via a new method, which comes from the maximum point of denominator of heat capacity. We also check the Maxwell equal area law in Sec. V. After that we calculate the critical exponents of the system in the extended phase space in Sec. VI. In Sec. VII, we examine the Ehrenfest equations at the critical point and confirm the validity of second order phase transition. In the last section we present our conclusions.

II. BASIC EQUATIONS
In the current paper, we set out to discuss the geometric and thermodynamic properties of charged black holes in d-dimensional GB-massive gravity with d − 4 compact dimensions. Regarding compactified extra dimensions, it has been shown that, depending on the horizon topology, one can obtain black string/membrane solutions in addition to black hole solutions. Furthermore, it has been pointed out that, in the context of GB gravity, one may obtain non-trivial modified solutions with an extra asymptotic charge [72]. In this paper, we focus on the black hole solutions with usual conserved charges.
The d-dimensional action of GB-massive gravity with the negative cosmological constant and in the presence of BI electrodynamics is where R, Λ, m, α and β are, respectively, the scalar curvature, the cosmological constant, the massive parameter, the GB factor and BI parameter. Also R µν and R µνγδ are Ricci and Riemann tensors, F = F µν F µν denotes the Maxwell invariant and f is a fixed symmetric tensor. In Eq. (1), c i 's are constants and U i 's are symmetric polynomials of the eigenvalues of d × d matrix K µ ν = √ g µα f αν , which can be written in the following forms Using the action (1) and variation of this action with respect to the metric tensor (g µν ) and the Faraday tensor (F µν ), respectively, lead to in the above equations G µν is the Einstein tensor, H µν and χ µν are and A. Black hole solutions Considering the metric of d-dimensional spacetime as where h ij dx i dx j is the line element with constant curvature (d − 2) (d − 3)κ and volume V d−2 and the ansatz metric in the following form [42] f µν = diag(0, 0, c 2 h ij ), (11) in which c is a positive constant. The metric function was obtained in Ref. [1] as where m 0 and q are integration constants which are related to the total mass and the electric charge of black hole, respectively. The notation d i is introduced by us to denote the term d−i (Recall that d is the spacetime dimensionality) so as to simplify the expressions of physical quantities in this paper. For example, d 1 denotes the term d − 1 while d 2 denotes the term d − 2. It is notable that, in the above solution, we used the gauge potential ansatz A µ = h(r)δ 0 µ in the Maxwell equation (7). Also, H, η and consistent h(r) are in the following forms It was shown that the asymptotical behavior of the solutions are (a)dS solutions with an effective cosmological constant (Λ ef f ) [1]. This effective cosmological constant reduces to ordinary Λ for vanishing α. It was also shown that neither massive nor BI parts affect the asymptotical behavior of the solutions [1].

B. Thermodynamics
The Hawking temperature of the black hole is given by [1] T = 1 4πN where N = 2ακd 3 d 4 + r 2 + and also, Υ ). It is notable that r + in the above expression denotes the largest real root of equation f (r) = 0. The total charge, the electric potential (U ) and the entropy of the black hole are [1] where H + = H r=r+ . Total mass of the black hole is in the following form [1] The first law of thermodynamics for black hole solution in the GB-BI-massive gravity was checked in Ref. [1] and it was found that these thermodynamics quantities satisfy the first law of black hole thermodynamics as

III. GEOMETRICAL THERMODYNAMICS
Here, we are interested in studying the critical behavior of the black holes through the use of geometrical method. This method builds phase space of black holes by using one of the thermodynamical quantities as thermodynamical potential and its corresponding extensive parameters as components of phase space. By doing so, a metric is obtained in which thermodynamical properties of the system are stored in its Ricci scalar. Divergencies of the Ricci scalar point out two important places in thermodynamical behavior of the system; whether system goes under second order phase transition or it meets a bound point. A bound point is where heat capacity/temperature meets a root. In other words, in bound points a limit for having physical system (positive temperature) is given. On the other hand, in phase transition point, heat capacity has a divergency, implying that there is a discontinuity in heat capacity. In place of this divergency, a second order phase transition takes place.
There are several methods for constructing phase space of black holes through thermodynamical quantities; Weinhold [7,8], Ruppeiner [9,10], Quevedo [11,12] and HPEM [13][14][15]. A successful method should cover all the mentioned points without any extra divergency for its Ricci scalar. Existence of extra divergency or mismatch between divergency of Ricci scalar and phase transition (or bound points) indicate that there is a case of anomaly. Recently, it was shown that employing Weinhold, Ruppeiner and Quevedo may lead to existence of anomaly [13][14][15]. To overcome the problems of other methods, HPEM metric was proposed. The structure of HPEM metric is where M X = ∂M/∂X and M XX = ∂ 2 M/∂X 2 . Now, by using total mass of black holes (20) with entropy (19) and electric charge (17), one can construct phase space and calculate its Ricci scalar. Due to economical reasons, we will not present obtained Ricci scalar but rather present its results in following diagrams (Figs. 1-3). Evidently, the number of phase transition points and their places are functions of massive ( Fig. 1), BI (Fig. 2) and GB (Fig. 3) parameters. For considered values of different parameters, these black holes enjoy the absence of bound point. In other words, for all values of the horizon radius, physical black holes exist. On other other hand, these black holes have second order phase transition in their thermodynamical structure. The number of these phase transitions may vary from one (see left panel of Fig. 1 and right panel of Fig. 3) to several (see Fig. 2) phase transitions.
The system has positive temperature but depending on the choices of different parameters, temperature may acquire one to several extrema. These extrema are where the heat capacity meets divergency. In other words, extrema of the temperature are places in which the heat capacity is divergent. Therefore, these extrema are places in which black holes go under the second order phase transition. The number of divergencies in the heat capacity, hence, phase transitions, is an increasing function of the massive (see Fig. 1) and BI (see Fig. 2) parameters while it is a decreasing function of the GB parameter (see Fig. 3).
Regarding BI theory, for large values of the nonlinearity parameter, system behaves like Reissner-Nordström black hole. Meaning for large values of this parameter, the effect of nonlinearity decreases and system behaves like it is in the presence of Maxwell theory of electromagnetic field. On the other hand, for small values of nonlinearity parameter, system has Schwarzschild like behavior. Taking these limits into account, one can extract following conclusions: The highest number of phase transitions, hence, the highest complexity in phase structure of these black holes is acquired for linear electromagnetic field. The generalization to nonlinear electromagnetic field reduces the number of phase transition and it may omit some of the phase transitions. By increasing power of the nonlinearity (decreasing the nonlinearity parameter), system would obtain the least number of phase transition which acquirable for charged black holes in this theory of the nonlinear electromagnetic field.
The GB gravity is a higher order gravity. In other words, value of the Ricci scalar, which is a parameter to measure curvature of the system, is higher in this theory of gravity comparing to Einstein gravity. Therefore, the gravity in this theory is stronger comparing to Einstein theory of gravity. The GB gravity provides an extra degree of freedom in term of GB parameter. Increasing GB parameter leads to increasing the value of Ricci scalar, hence, power of the gravity. We see for these black holes, that by increasing GB parameter, the number of phase transitions decreases. Meaning that for system with higher power of gravity (larger curvature), the number of phase transition and complexity in phase structure of these black holes decrease. Therefore, gravity here has an opposing effect on the number of phase transition.
Massive parameter is directly related to the mass of graviton. Plotted diagrams for variation of the massive parameter show that as the mass of graviton increases the black holes under consideration go under more phase transitions. In other words, by increasing the mass of graviton the complexity in thermal behavior and phase structure of these black holes increase.
It is evident that using HPEM metric provides suitable divergencies in its Ricci scalar for phase transitions that are observed in the heat capacity. In other words, divergences of the Ricci scalar of HPEM metric coincide with the phase transition points of heat capacity. Therefore, these two approaches yield consistent results. On the other hand, depending on the type of phase transition (smaller to larger or larger to smaller black holes), the sign of divergency of the Ricci scalar differs. If the transition is from larger to smaller, the sign of the Ricci scalar around the corresponding transition is positive, while the opposite (negative sign) is observed for the transition of smaller to larger black holes. These two differences in sign enable one to determine the type of phase transition of a system.

IV. P − V CRITICALITY THROUGH NEW APPROACH
In this section, we will regard critical behavior of these black holes through the use of a new method which was introduced in Ref. [73]. This method employs denominator of the heat capacity of black holes to extract a relation for pressure. This relation is independent from equation of state. The maximum of obtained relation is where phase transition takes place. In other words, the pressure and horizon radius of maximum of this relation is where system goes under the second order phase transition and Van der Waals like behavior is observed. In addition, the picture that this method draws for pressure smaller/larger than critical pressure is consistent with thermodynamical behavior for the system with same pressure in usual thermodynamical systems. In other words, for pressures smaller than critical pressure, two horizon radii exist, which marks two different phases in phase diagrams while for a pressure larger than critical pressure, no phase transition is observed. This is consistent with behavior of the T − V diagrams in which for pressures larger than critical pressure, no phase transition region exists. The consistency of this method with other methods was investigated in several papers [45,46,74,75]. Now, we will employ this method to obtain the critical pressure and horizon radius of these black holes. First, we use the proportionality between the cosmological constant and pressure with the heat capacity Using Eq. (16), we will obtain the denominator of heat capacity ∂T ∂S Q in the following form in which E is Now, by solving Eq. (25) with respect to pressure, a relation for pressure is obtained In order to study the critical behavior of these black holes, we should see whether a maximum exists for this relation. To do so, we employ numerical method. The results are presented in the following diagrams (Figs. 4 -7). First of all, it is evident that due to existence of maximum, these black holes enjoy a second order phase transition in their phase space. The critical pressure is a decreasing function of the massive, GB and BI parameters while their corresponding critical horizon radius are increasing functions of them (left panels of Figs. 4 -6). On the contrary, the critical horizon radius is a decreasing function of the dimension while the critical pressure is an increasing function of this parameter (left panel of Fig. 7).
Depending on choices of different parameters, one may come across two interestingly different behavior for P − r + diagrams; I) in one behavior, only one extremum exists for these diagrams (left panels of Figs. 4 and 7). II) in the other one, one minimum and one maximum exist (left panels of Figs. 5 and 6).
Considering mentioned concept for this method, only in maximum a second order phase transition exists. Therefore, we have a second order phase transition for both cases. On the other hand, for second behavior, for critical pressure, two horizon radii exist (due to formation of tail). This indicates that another branch for critical behavior exists. This critical behavior is not a second order phase transition but rather another kind (considering the concept of maximum). Interestingly, the second case of behavior for P − r + observed for small values of nonlinearity parameter and large values of GB parameter. In other words, for small values of BI parameter and large values of GB parameter, existence of extra branch in phase diagrams of these black holes is evident (left panels of Figs. 5 and 6). Existence of such behavior points out that another branch for phase diagrams exists for these black holes which is absent in other black holes. Such behavior is precisely due to existence of massive gravity. This means that by considering a massive theory of gravity for these black holes, another type of phase transition takes place. This emphasizes on the role and effects of massive gravity in thermodynamical behavior of these black holes.
In order to complete our study here, we will plot coexistence curves for variation of different parameters as well (right panels of Figs. 4 -7). The coexistence curves are representing small/larger black holes with similar pressure and temperature. The critical point is located at the end of this line which indicating after this point, phase transition does not take place. Evidently, the critical temperature is an increasing function of the massive gravity (right panel of Fig. 4) and dimensionality (right panel of Fig. 7) while it is a decreasing function of the GB (right panel of Fig.  5) and BI parameters (right panel of Fig. 6). Here, in these phase diagrams, we see that the presence of other phase transition is not observed. This indicates that our earlier interpretation is right. In other words, the branch for phase transition which was observed in P − r + diagram is not a second order phase transition. Also, we should point out that plotted diagrams indicate that no reentering of phase transition takes place for these black holes. The expressions of Hawking temperature and the entropy are listed in Eqs. (16) and (19) respectively. For T − r + graph, the possible critical point can be determined through ∂T ∂r + q=qc,r=r+c = 0, (28) For T − S graph, the possible critical point can be determined through Eqs. (30) and (31) are related to Eqs. (28) and (29) by Considering the above two relations and the fact that ∂S ∂r+ > 0, it is not difficult to conclude that the critical point conditions for T − r + and T − S graphs are equivalent to each other.
To probe the effect of massive gravity on critical quantities of T − S graph, we fix other parameters and let m vary from 0 to 0.3. The results are listed in Table I Table II. Lastly, we let β vary and keep other parameters fixed to study the effect of BI theory. The results are listed in Table III.     From Fig. 8, one can see clearly that both the graphs can be divided into three branches. The medium radius branch is unstable while both the large radius branch and the small radius branch are stable. The unstable branch in T − S curve can be removed with a bar vertical to the temperature axis T = T * as the approach in Ref. [76]. The possible Maxwell equal area laws for T − S and T − r + graphs read Note that S 1 , S 2 , S 3 denote the three values of entropy from small to large corresponding to T = T * while r 1 , r 2 , r 3 denote the three values of r + from small to large corresponding to T = T * .
To determine T * , one should first study the behavior of free energy (F ), which can be obtained as We plot the free energy for the case q = 0.2q c , m = 0.2, α = 0.5, β = 0.5, c = c 1 = c 2 = 2, c 3 = 0.2, c 4 = −0.2, d = 6, Λ = −0.1 in Fig. 9, where we can find the swallow tail characteristic of first order phase transition. Numerical  Table. IV and V, the relative errors are very small and the Maxwell equal area law holds for not only T − r + curves, but also T − S curves.     The possible Maxwell equal area laws for P − r + and P − V graphs read Here, r 1 , r 2 and r 3 denote the three values of r + from small to large corresponding to P = P * in P − r + graph while V 1 , V 2 and V 3 denote the three values of V from small to large corresponding to P = P * in P − V graph. Note that thermodynamic volume V is defined in the extended phase space as V = We plot the Gibbs free energy for the case P * = 0.5P c , m = 0.5, Fig. 10. The classical swallow tail behavior can also be found. We further calculate both the left hand side and right hand side of Eqs. (37) and (38). As shown in Tables. VI and VII, the relative errors for P − r + graph are very large while those for P − V graph are amazingly small, leading to the conclusion that the Maxwell equal area law holds for P − V graph while it fails for P − r + graph. Our numerical results here for the GB-BI-massive black holes further backup the findings in former researches [77,78].

VI. CRITICAL EXPONENTS
To characterize the critical behavior near the critical point in the extended phase space, one usually introduces the following critical exponents Note that we use the notations α 1 and β 1 instead of the classical notation α and β here because α and β already have other meanings in this paper. As can be seen from the above definitions, the exponents α 1 , β 1 , γ and δ describe the behavior of specific heat C V , the order parameter η, the isothermal compressibility coefficient κ T and the critical isotherm respectively.
Before calculating the above critical exponents, it would be convenient to define The equation of state in the extended phase space has been derived in Ref. [1] as where . Identifying the specific volume v as v = 4r+ d2 , the equation of state can be reorganized as Then the equation of state in the extended phase space can be expanded as p = 1 + p 10 t + p 01 ǫ + p 11 tǫ + p 02 ǫ 2 + p 03 ǫ 3 + O(tǫ 2 , ǫ 4 ).
where the expansion coefficients can be calculated as From the equal area law, one can further derive ǫs ǫ l ǫ dp dǫ dǫ = 0, where dp dǫ can be calculated as p 11 t + 3p 03 ǫ 2 . Denoting the subscript "l" and "s" as the quantity of large black hole and small black hole, respectively, one can obtain On the other hand, the pressure of large black hole equals to that of small black hole as follow 1 + p 10 t + p 11 tǫ l + p 03 ǫ 3 l = 1 + p 10 t + p 11 tǫ s + p 03 ǫ 3 s , because during the phase transition the pressure of the black hole keeps unchanged. With Eqs. (53) and (54), one can get So the order parameter can be derived as leading to the conclusion that β 1 = 1/2. It is not difficult to deduce that with which one can draw the conclusion that γ = 1.
One can obtain the critical isotherm by substituting t = 0 into Eq. (47) implying that δ = 3. The entropy S does not depend on the Hawking temperature T . So the specific heat with fixed volume C V is equal to zero, with the critical exponent α 1 = 0.
The above exponents are totally the same as those in former literature. It can be attributed to the effect of mean field theory.

VII. ANALYTICAL CHECK OF EHRENFEST EQUATIONS AT THE CRITICAL POINT IN THE EXTENDED PHASE SPACE
It is important to classify the nature of phase transition. As we know, the Clausius-Clapeyron equation is satisfied for a first order phase transition while for a second order phase transition one can utilize the famous Ehrenfest equations as follows where the volume expansion coefficientα = 1 V ( ∂V ∂T ) P and isothermal compressibility coefficient κ T = − 1 V ( ∂V ∂P ) T . Note that we use the notationα instead of the classical notation α here because α already has other meaning in this paper.
Utilizing the definition ofα, one can derive So the R.H.S of Eq. (59) can be obtained as The subscript "c" here denotes the corresponding quantity at the critical point. It is not difficult to obtain ∆C P T V ∆α = d 2 (r 2 c + 2d 4 d 3 ακ) From Eqs. (63) and (64), we can draw the conclusion that the first equation of Ehrenfest equations is valid at the critical point.
The L.H.S of Eq. (60) can be obtained as With both the definitions of κ T andα, one can deduce The above equation and the validity of Ehrenfest equations show that GB-BI-massive black holes undergo second order phase transition at the critical point of P − V criticality in the extended phase space. The result here is consistent with the nature of liquid-gas phase transition at the critical point and support the findings in former literatures [79][80][81].

VIII. CLOSING REMARKS
In this paper, we have studied thermodynamical behavior of Einstein-GB-massive black holes in the presence of BI nonlinear electromagnetic field near critical point.
First, some comments regarding the effects of mass of graviton, nonlinearity of the electromagnetic field and power of gravity (value of the GB curvature term) on phase structure and its complexity were given. In addition, geometrical thermodynamics was used to investigate phase transition of these black holes based on canonical ensemble.
Next, by using the denominator of heat capacity and the proportionality between the cosmological constant and thermodynamic pressure, critical behavior of these black holes was investigated. It was shown that these black holes enjoy an anomaly in their phase structure. In other words, in addition to Van der Waals like phase transition in their phase diagrams, these black holes enjoy another type of phase transition which is different from usual Van der Waals like phase transition. Plotted coexistence curves also confirmed that only one second order phase transition exists for these black holes.
Moreover, Maxwell equal area law was employed to investigate the Van der Waals like behavior and structure of these black holes. It was shown that Maxwell equal area law holds for T − r + , T − S and P − V diagrams while it fails regarding P − r + curves. Calculations regarding critical exponent showed that these exponents are independent of massive gravity and are same as those derived previously. Finally, the Ehrenfest equations were used to determine the type of phase transition. It was shown that these black holes undergo second order phase transition at the critical point.