Einstein-Born-Infeld-Massive Gravity: adS-Black Hole Solutions and their Thermodynamical properties

In this paper, we study massive gravity in the presence of Born-Infeld nonlinear electrodynamics. First, we obtain metric function related to this gravity and investigate the geometry of the solutions and find that there is an essential singularity at the origin ($r=0$). It will be shown that due to contribution of the massive part, the number, types and places of horizons may be changed. Next, we calculate the conserved and thermodynamic quantities and check the validation of the first law of thermodynamics. We also investigate thermal stability of these black holes in context of canonical ensemble. It will be shown that number, type and place of phase transition points are functions of different parameters which lead to dependency of stability conditions to these parameters. Also, it will be shown how the behavior of temperature is modified due to extension of massive gravity and strong nonlinearity parameter. Next, critical behavior of the system in extended phase space by considering cosmological constant as pressure is investigated. A study regarding neutral Einstein-massive gravity in context of extended phase space is done. Geometrical approach is employed to study the thermodynamical behavior of the system in context of heat capacity and extended phase space. It will be shown that GTs, heat capacity and extended phase space have consistent results. Finally, critical behavior of the system is investigated through use of another method. It will be pointed out that the results of this method is in agreement with other methods and follow the concepts of ordinary thermodynamics.

In this paper, we study massive gravity in the presence of Born-Infeld nonlinear electrodynamics. First, we obtain metric function related to this gravity and investigate the geometry of the solutions and find that there is an essential singularity at the origin (r = 0). It will be shown that due to contribution of the massive part, the number, types and places of horizons may be changed. Next, we calculate the conserved and thermodynamic quantities and check the validation of the first law of thermodynamics. We also investigate thermal stability of these black holes in context of canonical ensemble. It will be shown that number, type and place of phase transition points are functions of different parameters which lead to dependency of stability conditions to these parameters. Also, it will be shown how the behavior of temperature is modified due to extension of massive gravity and strong nonlinearity parameter. Next, critical behavior of the system in extended phase space by considering cosmological constant as pressure is investigated. A study regarding neutral Einsteinmassive gravity in context of extended phase space is done. Geometrical approach is employed to study the thermodynamical behavior of the system in context of heat capacity and extended phase space. It will be shown that GTs, heat capacity and extended phase space have consistent results. Finally, critical behavior of the system is investigated through use of another method. It will be pointed out that the results of this method is in agreement with other methods and follow the concepts of ordinary thermodynamics.

I. INTRODUCTION
Regarding experimental agreements of Einstein gravity (EN) in various area of astrophysics and cosmology, motivates one to consider it as an acceptable theory. In addition, adding a constant term Λ in the EN-Hilbert action may lead to agreement between the results of EN-Λ gravity with dark energy prediction.
On the other hand, general relativity is consistent with interaction of massless spin 2 fields, in which related gravitons are massless particles with two degrees of freedom. Since the quantum theory of massless gravitons is non-renormalizable [1], one may be motivated for modifying general relativity to massive gravity. In order to build up a massive theory with a massive spin 2 particle propagation, one can add a mass term to the EN-Hilbert action. This will result into graviton having a mass of m which in case of m → 0, the effect of massive gravity will be vanished. A class of massive gravity theory in flat and curved background which leads to absence [2] and existence [3] of ghost, have been investigated. Also, the quantum aspects of the massive gravity and a nonlinear class of massive gravity in ghost-free field [4] have been explored in Refs. [5,6]. Generalization to nonlinearly charged massive black holes was done in Ref. [7]. More details regarding the motivations of massive gravity is given in Ref. [8].
In this paper, we are interested in studying the nontrivial adS massive theory that was investigated in [9,10]. The motivation for this consideration is due to fact that graviton shows similar behavior as lattice in holographic conductor [11]. In other words, a Drude like behavior is observed for the case of massless graviton in this theory which makes the role of graviton similar to lattice. Another interesting subject for study in this theory is metal-insulator transition [12]. Recently, charged massive black holes with consideration of this theory have been investigated in [13]. The P − V criticality of these solutions and their geometrical thermodynamic aspects have been studied [14,15]. Also, the generalization to Gauss-Bonnet-Maxwell-massive gravity and its stability, geometrical thermodynamics and P − V criticality have been investigated [16].
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 hence, it has infinite self-energy. To overcome this problem in classical electrodynamics, Born and Infeld in Ref. [17] introduced a nonlinear electromagnetic field, with main motivation, to solve infinite self-energy problem by imposing a maximum strength of the electromagnetic field. Then, Hoffmann in Ref. [18] investigated EN gravity in the presence of Born-Infeld (BI) electrodynamics. In recent two decades, exact solutions of gravitating black objects in the presence of BI theory have been vastly investigated [19,20]. Another interesting property of BI is that, BI type effective action arises in an open superstring theory and D-branes are free of physical singularities [21]. For a review of aspects of BI theory in the context of string theory see Ref. [22]. Recently, there has been growing interest in Eddington-inspired BI gravity in context of black holes and cosmology [23]. Also, it was proposed that one can consider BI theory as a gravitational theory [24]. Dualization of the BI theory and some of the special properties of this theory have been investigated in Ref. [25].
There are several approaches for studying and obtaining critical behavior and phase transition points of black holes: First method is based on studying heat capacity. It was pointed out that roots and divergencies of the heat capacity are representing phase transition points. In other words, in place of roots and divergencies of the heat capacity system may go under phase transition. Another important property of the heat capacity is investigation of the thermal stability. Systems with positive heat capacity are denoted to be in thermally stable states. Therefore, the stability conditions are indicated by changes in sign of heat capacity [26]. This is known as canonical ensemble.
In the second method, by using the renewed interpretation of cosmological constant as thermodynamical variable, one can modify the thermodynamical structure of the phase space [27]. One of the most important property of this method is the similarity of critical behavior of the black holes and ordinary thermodynamical Van der Waals liquid/gas systems [28]. Recently, it was pointed out that the extended phase space should be interpreted as an RG-flow in the space of field theories, where isotherm curves codify how the number of degrees of freedom N (or the central charge c) runs with the energy scale [29]. On the other hand, it was shown that variation of cosmological constant could be corresponded to variation of number of the colors in Yang-Mills theory residing on the boundary spacetime [30].
The third method is using geometrical concept for studying critical behavior. In other words, by employing a thermodynamical potential and its corresponding extensive parameters, one can construct phase space. The divergencies of Ricci scalar in constructed metric are denoted as phase transition points. There are several metrics for this method that one can name: Weinhold [31], Quevedo [32] and HPEM [33] which has mass as thermodynamical potential and Ruppeiner [34] in which entropy is considered as thermodynamical potential. These metrics are used in context of heat capacity. Another set of metrics was introduced in Ref. [35] which can be used in context of extended phase space.
Finally, a fourth method was introduced in Ref. [35] which is based on denominator of the heat capacity. In this method by replacing cosmological constant with pressure in denominator of the heat capacity and solving it with respect to pressure, a new relation is obtained for pressure. The existence of maximum for obtained relation, represents the critical pressure and volume in which phase transition takes place. The behavior of system in case of this method is consistent with ordinary thermodynamical concepts [15,35].
The outline of the paper will be as follow. In Sec. II, we introduce action and basic equations related to EN-BI-massive gravity. Sec. III is devoted to obtain the black hole solutions of this gravity and investigation of the geometrical structure of them. In the next section, we calculate conserved and thermodynamic quantities related to obtained solutions and check the validation of the first law of thermodynamics. In section V, we study thermal stability of the EN-BI-massive black hole solutions in canonical ensemble. Next, we consider cosmological constant as pressure and study the critical behavior of the system. Then we employ the geometrical methods for investigating thermodynamical behavior of the system and extend this study by another method. In the last section we present our conclusions.

II. BASIC EQUATIONS
The d-dimensional action of EN-massive gravity with negative cosmological constant and a nonlinear electrodynamics is where R is the scalar curvature, is the negative cosmological constant and f is a fixed symmetric tensor. In Eq. (1), c i are constants and U i are symmetric polynomials of the eigenvalues of the d × d matrix K µ ν = √ g µα f αν which can be written as follows Here, we want to study a particular model of nonlinear electrodynamics called BI theory which has attracted lots of attentions due to its relation to effective string actions. The function L(F ) for BI theory is given as where β is the BI parameter, the Maxwell invariant is F = F µν F µν in which F µν = ∂ µ A ν −∂ ν A µ is the electromagnetic field tensor and A µ is the gauge potential. Variation of the action (1) with respect to the metric tensor g µν and the Faraday tensor F µν , leads to where G µν is the EN tensor and χ µν is the massive term with the following form

III. BLACK HOLE SOLUTIONS IN EN-BI-MASSIVE GRAVITY
In this section, we obtain static nonlinearly charged black holes in context of massive gravity with adS asymptotes. For this purpose we consider the metric of d-dimensional spacetime in the following form where h ij dx i dx j is a (d − 2) dimension line element for an Euclidian space with constant curvature (d − 2) (d − 3)k and volume V d−2 . We should note that the constant k, which indicates that the boundary of t = constant and r = constant, can be a negative (hyperbolic), zero (flat) or positive (elliptic) constant curvature hypersurface. We consider the ansatz metric [13] f µν = diag(0, 0, c 2 h ij ), where c is a positive constant. Using the metric ansatz (7), U i 's are in the following forms [13] Using the gauge potential ansatz A µ = h(r)δ 0 µ in electromagnetic equation (4) and considering the metric (6), we obtain in which H is the following hypergeometric function where Γ = d2d3q 2 β 2 r 2d 2 and q is an integration constant which is related to the electric charge. Also, the electromagnetic field tensor in d-dimensions is given by Now, we are interested in obtaining the static black hole solutions. One may use components of Eq. (3) and obtain metric function f (r). We use the tt and x 1 x 1 components of the Eq. (3), which can be written as We can obtain the metric function f (r), by using the Eqs. (11) and (12) with the following form where m 0 is an integration constant which is related to the total mass of the black hole. It should be noted that, obtained metric function (13), satisfy all components of the Eq. (3), simultaneously. Now, we are in a position to review the geometrical structure of this solution, briefly. We first look for the essential singularity(ies). The Ricci scalar and the Kretschmann scalar are lim and so confirm that there is a curvature singularity at r = 0. The Ricci and Kretschmann scalars are 2d d2 Λ and 8d d1d 2 2 Λ 2 at r −→ ∞. Therefore, the asymptotic behavior of these solutions are (a)dS for (Λ < 0 ).
On the other hand, in the absence of massive parameter (m = 0), the solution (13) reduces to an d-dimensional asymptotically adS topological black hole with a negative, zero or positive constant curvature hypersurface in the following form In order to study the effects of the EN-BI-massive gravity on metric function, we have plotted various diagrams (Figs. 1 -3 ).
By considering specific values for the parameters, metric function has different behaviors. Depending on the choices of the parameters, EN-BI-massive black holes can behave like Reissner-Nordström black holes. In other words, these black holes may have two horizons, one extreme horizon and without horizon (naked singularity) (see Fig. 1 for more details). On the other hand, by adjusting some of the parameters of EN-BI-massive black holes, we encounter with interesting behaviors. The solutions may have three or higher horizons (Figs. 2 and 3). The existence of three or higher horizons for black holes is due to the presence of massive gravity [15,16]. In addition to the significant effects of massive term, we should note that the nonlinearity parameter can affect on the number of horizons. In addition, β can change the type of singularity. In other words, depending on the parameters, one can find a β c in which singularity is spacelike for β < β c , and it would be timelike for β > β c (see [20] for more details). Now, we give a brief discussion regarding Carter-Penrose diagrams. In order to study the conformal structure of the solutions, one may use the conformal compactification method through plotting the Carter-Penrose (conformal) diagrams. As we mentioned before, depending on the value of nonlinearity parameter, β, one may encounter with timelike or spacelike singularity. Penrose diagrams regarding to timelike singularity was discussed in [16] (see conformal diagrams in [16]). Here we focus on special case in which singularity is spacelike (β < β c ). In other words, the singularity of this nonlinearly charged black holes behaves like uncharged Schwarzschild solutions (see Fig. 3). This means that, although massive and nonlinearity parts of the metric function can change the type of singularity and horizon structure of black holes, they does not affect asymptotical behavior of the solutions. Drawing the Carter-Penrose diagrams, we find the causal structure of the solutions are asymptotically well behaved.

IV. THERMODYNAMICS
In this section, we calculate the conserved and thermodynamic quantities of the static EN-BI-massive black hole solutions in d-dimensions and then check the first law of thermodynamics.
By using the definition of Hawking temperature which is related to the definition of surface gravity on the outer horizon r + , one can find where Γ + = d2d3q 2 β 2 r 2d 2 + . The electric charge, Q, can be found by calculating the flux of the electric field at infinity, yielding In order to obtain the entropy of the black holes, one can employ the area law of the black holes. It is a matter of calculation to show that entropy has the following form [36] It was shown that by using the Hamiltonian approach or counterterm method, one can find the mass M of the black hole for massive gravity as in which by evaluating metric function on the horizon (f (r = r + ) = 0), one can obtain where , −Γ + . It is notable that, U is the electric potential, which is defined in the following form Having conserved and thermodynamic quantities at hand, we are in a position to check the first law of thermodynamics for our solutions. It is easy to show that by using thermodynamic quantities such as charge (18), entropy (19) and mass (20), with the first law of black hole thermodynamics we define the intensive parameters conjugate to S and Q. These quantities are the temperature and the electric potential which are the same as those calculated for temperature (17) and electric potential (22).

V. HEAT CAPACITY AND STABILITY IN CANONICAL ENSEMBLE
Here, we study the stability conditions and the effects of different factors on these conditions. The stability conditions in canonical ensemble are based on the signature of the heat capacity. The negativity of heat capacity represents unstable solutions which may lead to following results: unstable solutions may go under phase transition and acquire stable states. This phase transition could happen whether when heat capacity meets a root(s) or has a divergency. Therefore, the roots of regular numerator and denominator of the heat capacity are phase transition points. In the other scenario, the heat capacity is always negative. This is known as non-physical case. But there is a stronger condition which is originated from the temperature. The positivity of the temperature represents physical solutions whereas its negativity is denoted as non-physical one. Therefore, in order to getting better picture and enriching the results of our study, we investigate both temperature and heat capacity, simultaneously.
The heat capacity is given by Considering Eqs. (17) and (19), it is a matter of calculation to show that In order to study the effects of different parameters on stability conditions and temperature, we have plotted various diagrams .
Interestingly, in the absence of the massive parameter ( Fig. 4 right panel), temperature starts from −∞ and it is only an increasing function of horizon radius with a root. Therefore, we have two regions of physical and non-physical solutions. Adding massive gravity could modify the behavior of the temperature into an U shape diagram starting from +∞ without any root. The extremum is an increasing function of massive parameter (Fig. 4 right panel) and dimensions ( Fig. 8 right panel), whereas, it is a decreasing function of k (Fig. 6 right panel) and electric charge ( Fig. 7 right panel). The only exception for this behavior is for strong nonlinearity parameter. Interestingly, for large values of nonlinearity parameter, a massive-less like behavior is observed (Fig. 5 right panel). In other words, the temperature starts from −∞ but the effect of the massive could be seen through two extrema. The root of temperature and smaller extremum are increasing functions of β and larger extremum is a decreasing function of it. Another interesting property of temperature is the effect of dimensions. Studying Fig. 8 (right panel) shows that the temperature for every two sets of dimensions will coincide with each other. In other words, there are places in which despite differences in dimensions, two black holes with two different dimensions and same values for other parameters will have same temperature in a special r + .
In absence of massive gravity, black holes could acquire temperature from zero to +∞, whereas adding massive, will cause the black holes never acquire some temperature. This effect is vanished in case of large nonlinearity parameter. In other words, the strength of nonlinearity parameter has opposing effects to massive's ones. Also, the U shape diagram indicates that for every temperature that black holes can acquire two horizons exist except for the extremum. Therefore, considering Hawking radiation, one is not able to recognize the size of these black holes by measuring their Hawking radiation. It is worthwhile to mention that extrema and root(s) of temperature are phase transition points of heat capacity.
Regarding stability, it is evident that in absence of the massive gravity, there exists a region of the instability which is located where the temperature is negative. Therefore, this is a non-physical solution (Fig. 4 left panel). Interestingly, by adding massive gravity, the non-physical region is vanished and heat capacity acquires divergence point without any root. Before divergence point, the heat capacity is negative. Therefore, in this region black holes are unstable. In divergence point, black holes go under phase transition of smaller unstable to larger stable black holes. The divergence point is an increasing function of massive parameter (Fig. 4 left panel) and dimensions (Fig. 8  left panel), whereas, it is a decreasing function of k (Fig. 6 left panel) and electric charge (Fig. 7 left panel). Interestingly, in strong nonlinearity parameter, the mentioned behavior is modified. In this case black holes enjoy one root and two divergence points. Before root and between two divergencies, heat capacity is negative and between root and smaller divergence point and after larger divergence point, heat capacity is positive. According to thermodynamical concept, systems go under phase transition to acquire stable states. Therefore, following phase transitions take place: non-physical unstable to physical stable (in root), large unstable to smaller stable (in smaller divergence point) and smaller unstable to larger stable black holes (in larger divergency). Root and smaller divergence point are increasing functions of β (Fig. 5 left panel), whereas, larger divergency is a decreasing function of it (Fig. 5 middle panel). It is worthwhile to mention that larger divergency is not highly sensitive to variation of nonlinearity parameter.
Comparing obtained results for heat capacity (regarding phase transitions) and the behavior of the temperature, one can see that larger to smaller phase transition takes place at maximum (compare Fig. 6 left diagram with right) and smaller to larger one happens at minimum (compare Fig. 6 middle diagram with right) of temperature. Therefore, one is able to recognize the type and number of phase transition by only studying temperature's diagrams.

VI. P − V CRITICALITY OF CHARGED BLACK HOLES IN EN-BI-MASSIVE GRAVITY
Now, we are in a position to study the critical behavior of the system through phase diagrams. Using the renewed interpretation of the cosmological constant as thermodynamical pressure, one can use following relation to rewrite thermodynamical relations of the solutions in spherical horizon [28] which results into following conjugating thermodynamical variable corresponding to pressure [28] V = ∂H ∂P S,Q .
Due to existence of the pressure in obtained relation for total mass of the black holes, one can interpret the total mass as thermodynamical quantity known as Enthalpy. This interpretation will lead to the following relation for Gibbs free energy [28] Now by using Eqs. (20) and (27) with the relations of volume and Gibbs free energy (Eqs. (28) and (29)), one finds and Obtained relation for volume indicates that volume of the black holes is only a function of the topology of the solutions and independent of electrodynamics and gravitational extensions, directly.
In order to obtain critical values, one can use P −V diagrams. In other words, by studying inflection point properties one can obtain critical values in which phase transitions may take place. Therefore, one can use Considering obtained values for temperature (17) and pressure (27), one can obtain pressure as Now, by considering Eq. (32) with obtained relation for pressure (33), one can obtain two relations for finding critical quantities. Due to economical reasons, we will not present them. Regarding the contribution of electromagnetic part, it is not possible to obtain critical horizon analytically, and therefore, we use numerical method. Considering the variation of β and massive parameter, one can draw following tables In addition, we plot following diagrams  to investigate that obtained values are the ones in which phase transition takes place or not.
The formation of swallow tail in G − T diagrams for pressure smaller than critical pressure (Figs. 9 and 11 right panels), subcritical isobars in T − r + diagrams for critical pressure (Figs. 9 and 11 middle panel) and isothermal diagrams in case of critical temperature in P − r + diagrams (Figs. 9 and 11 left panels), show that obtained values are critical ones in which phase transition takes place. It is evident that critical pressure (Fig. 10 left panel) and temperature (Fig. 10 middle panel) are increasing functions of the massive parameter, whereas the critical horizon ( Fig. 10 left and middle panels) and universal ratio of Pcrc Tc are decreasing functions of this parameter. It is worthwhile to mention that length of subcritical isobars (which is known as phase transition region) is a decreasing function of massive parameter (Fig. 10 middle panel). In opposite, the size of swallow tail and the energy of different phases are increasing functions of m (Fig. 10 right panel).
Interestingly, the effects of variation of nonlinearity parameter is opposite of massive parameter. In other words, critical pressure (Fig. 12 left panel), temperature (Fig. 12 middle panel) and the size of swallow tail (Fig. 12 left panel) are decreasing functions of β, whereas, the critical horizon radius (Fig. 12 left and middle panels), length of subcritical isobars (Fig. 12 middle panel) and universal ration of Pcrc Tc are increasing functions of nonlinearity parameter.
It should be pointed that the length of subcritical isobars affects single regions of different states which in our cases are smaller and larger black holes. In other words, increasing the length of subcritical isobars (phase transition region) decreases the single state regions.

A. Neutral Massive black holes
In this section, by cancelling the electric charge (q = 0), we will study the critical behavior of the system. Previously, it was shown that Schwarzschild black holes does not have any phase transition in context of extended phase space. Now, we are investigating the effects of massive gravity in case of EN-massive gravity. Using obtained relation for calculating critical horizon radius in previous part and setting q = 0, one can find following relation for calculating critical horizon radius It is a matter of calculation to show that this relation has following roots which are critical horizon radii Obtained relation shows that in absence of massive gravity, critical horizon radius will be zero which is not of our interest. This result consistent with Schwarzschild case. Now, for the simplicity, we consider the case of c 4 = 0. This leads into following critical horizon radius It is evident that for the cases of d = 4 and d > 4 with vanishing c 3 , the critical horizon radius will be zero. Therefore, there is no phase transition for these black holes. Interestingly for case of d > 4, the condition for having a positive critical horizon radius will be c 3 < 0 and 1 + c 2 c 2 m 2 > 0. By employing obtained value for critical horizon radius, one can find critical temperature and pressure in the following forms Considering obtained values, one can show that following equality is hold which shows that in this case, Pccrcc Tcc is a function of massive parameter and coefficients. Using obtained critical values (Eqs. (36), (37) and (38)) with Eqs. (17), (27), (29), (33) and setting q = 0, we plot following diagrams for 5 and 6 dimensions (Figs. 13 and 14).
In Ref. [37], it was shown that in context of neutral Gauss-Bonnet black holes, no phase transition is observed in 6-dimensions. Here, the extension of the massive gravity enables the black holes to enjoy the existence of second order phase transition in 6-dimensions. Also, we observed that contrary to Gauss-Bonnet case, Pccrcc Tcc is a function of massive gravity.

VII. GEOMETRICAL PHASE TRANSITION IN CONTEXT OF HEAT CAPACITY AND EXTENDED PHASE SPACE
In this section, we employ the geometrical concept for studying thermodynamical behavior of the obtained solutions. In order to do so, we employ HPEM method. In this method, the thermodynamical phase space is constructed by considering mass of the black holes as thermodynamical potential. By doing so, the components of the phase space will be extensive parameters such as electric charge, entropy and etc. The general form of HPEM metric is [33] where M S = ∂M/∂S, M SS = ∂ 2 M/∂S 2 and χ i (χ i = S) are extensive parameters which are components of phase space. Now, we will investigate whether the phase transition points that were obtained in section (IV ) coincide with all divergencies of the Ricci scalar of HPEM metric. For economical reasons, we only plot diagrams correspond to variation of massive and nonlinearity parameters. To do so, we use Eqs. (18), (19) and (20) with HPEM metric (Eq. 40). This leads into following diagrams (Fig. 15). It is evident that employed metric has consistent results with what were found in case of heat capacity (compare Figs. 4 and 5 with 15). In other words, the divergencies of the Ricci scalar are matched with phase transition points of the heat capacity. An interesting characteristic behavior of the diagrams is the different divergencies for different types of phase transition. In case of larger to smaller black holes phase transition, the divergency of the Ricci scalar is toward +∞ (compare behavior enables us to recognize the type of phase transition independent of heat capacity.
Next, we employ another geometrical metric for studying the critical behavior of the system in context of extended phase space. In this metric, Due to consideration of the cosmological constant as thermodynamical pressure, we have three extensive parameters; electric charge, entropy and pressure. In order to construct phase space we employ following metric [35] Considering Eqs. (20), (26), (27) and (41), we plot following diagram ( Fig. 16) with respect to Fig. 9. Due to existence of a root for heat capacity, in all plotted diagrams, a divergency is observed (Fig. 16 left panel). It is evident that for pressures smaller than critical pressure, system goes under two phase transitions with different horizon radii (Fig. 16 middle panel). This is consistent with what was observed in studying T − r + diagrams of Fig. 9. On the other hand, for critical pressure system goes under a phase transition. The place of this divergency is exactly located at the critical horizon which is obtainable through T − r + diagrams of Fig. 9 (Fig. 16 right panel). Finally, for pressures larger than critical pressure no phase transition is observed and the behavior of Ricci scalar will be what is plotted in Fig. 16 (left panel). These results are consistent with ordinary thermodynamical concepts and indicates that these three pictures (phase diagrams, heat capacity and geometrical thermodynamics) are in agreement.

VIII. HEAT CAPACITY AND CRITICAL VALUES IN THE EXTENDED PHASE SPACE
The final section of this paper is devoted to calculation of the critical pressure in extended phase space by using denominator of the heat capacity. It was shown that one can calculate critical pressures that were obtained in section (V ) by using denominator of the heat capacity [35]. To do so, one should replace the cosmological constant in denominator of the heat capacity (26) with its corresponding pressure (27). Then, solve the denominator of the heat capacity with respect to pressure. This will lead into following relation Obtained relation for pressure is different from what was obtained through use of temperature (33). In this relation, the maximum(s) of pressure and its corresponding horizon radius are critical pressure and horizon radius in which phase transition takes place. Now, by using indicated values in table 1 and Eq. (42), we plot following diagram (Fig. 17).
It is evident that obtained maximums are critical pressures in which phase transitions take place. The thermodynamical concept that was mentioned in last section (pressure being smaller than critical pressure leads to existence of two phase transitions and for pressures larger than critical pressure no divergency is observed) is also hold in case of this approach. In other words, this approach is an additional method for studying critical behavior of the system and the results of this approach is consistent with GTs, heat capacity and extended phase space.

IX. CONCLUSIONS
In this paper, we have considered EN-massive gravity in presence of BI nonlinear electromagnetic field. It was shown that considering this configuration leads to modification of the number and place of horizons that black holes can acquire. In other words, cases of multiple horizons were observed with different phenomenologies. Next, conserved and thermodynamical quantities were obtained and it was shown that first law of thermodynamics hold for these black holes.
Next, we studied the thermodynamical behavior of the system. It was shown that temperature in the presence and absence of massive gravity presents different behaviors. Adding massive put limitations on values that temperature can acquire, while, there was no limitation for temperature in the absence of it. Interestingly, this behavior was modified in the presence of strong nonlinearity parameter. In strong nonlinearity parameter, the behavior of temperature returned to a massive-less like behavior but the effects of the massive were observed in existences of extrema. It was also seen, that in case of different dimensions, for each pair of dimensions, one can find a point in which temperature for both dimensions are equal.
Regarding the stability, it was seen that in the presence of massive gravity, black holes enjoy one phase transition of the smaller unstable to larger stable. The phase transition was related to the divergency of heat capacity. Then again, in strong nonlinearity parameter, this behavior was modified. In this case, black holes had three phase transitions of smaller non-physical unstable to larger physical stable (in place of root), larger unstable to smaller stable (in place of smaller divergency) and smaller unstable to larger stable (in place of larger divergency).
Clearly, one can conclude that nonlinear electromagnetic field has an opposing effect comparing to massive gravity. Strong nonlinearity parameter modifies the effects of the massive gravity and return the system to the massive-less like behavior, although the effects of massive still observed through extrema.
It was pointed out that at maximums of the temperature, larger unstable to smaller stable and at minimums, smaller unstable to larger stable phase transitions take place. Therefore, studying temperature provides an independent picture for studying phase transitions and stability of the solutions.
Next, we extended phase space by considering cosmological constant as thermodynamical variable known as pressure. It was shown that volume of the black holes is independent of generalization of the electromagnetic field and extension of the massive gravity. Obtained values where critical points in which phase transitions took place. It was shown that the effects of variation of nonlinearity parameter was opposite of the massive parameter. In other words, these two factors put restrictions on each others effects.
Interestingly, in Ref. [15] variation of massive gravity highly modified the critical temperature and pressure. In case of obtained solutions in this paper, the modification was not as considerable as what was observed in case of Gauss-Bonnet-Maxwell-massive black holes. This shows that generalization of electromagnetic field puts stronger restrictions on the effects of the massive gravity. In other words, in order to have stronger control over contributions of the massive gravity one should increase the nonlinearity of the electromagnetic sector.
In addition, a study in context of neutral solutions was conducted. It was shown that due to contribution of the massive gravity, the chargeless solutions of this gravity also enjoy the existence of phase transition. In other words, black holes in EN-massive gravity go under phase transitions in extended phase space. Also, it was shown that ratio of Pcrc Tc was a function of massive gravity. It was also pointed out that in 6-dimensions, contrary to case of Gauss-Bonnet black holes, these black holes enjoy second order phase transition. In addition, it was shown that in case of d = 4, no phase transition for massive black holes is observed.
Next, geometrical approach was used for studying critical behavior of the system in context of heat capacity and extended phase space. It was shown that employed metrics for both cases have consistent results and follow the concepts of ordinary thermodynamics. The characteristic behavior of divergencies in Ricci scalar of the geometrical thermodynamical metrics, enabled us to recognize the type of phase transition (smaller to larger or larger to smaller).
Finally, another method which was based on denominator of the heat capacity was used to calculate critical pressure and horizon radius. It was shown that this method has consistent results with extended phase space and follow the concepts of ordinary thermodynamics. In other words, this method provides an independent approach for investigating critical behavior of the system. Due to generalization of Born-Infeld for electromagnetic sector of the solutions, it will be worthwhile to study the effects of this generalization on conductivity of these black holes and their corresponding superconductors phase transition. Specially, it will be interesting to see how this generalization will affect the interpretation of graviton as lattice and the Drude like behavior. Also, it will be worthwhile to study the metal-insulator transition in context of these solutions.