Massive charged BTZ black holes in asymptotically (a)dS spacetimes

Motivated by recent developments of BTZ black holes and interesting results of massive gravity, we investigate massive BTZ black holes in the presence of Maxwell and Born-Infeld (BI) electrodynamics. We study geometrical properties such as type of singularity and asymptotical behavior as well as thermodynamic structure of the solutions through canonical ensemble. We show that despite the existence of massive term, obtained solutions are asymptotically (a)dS and have a curvature singularity at the origin. Then, we regard varying cosmological constant and examine the Van der Waals like behavior of the solutions in extended phase space. In addition, we employ geometrical thermodynamic approaches and show that using Weinhold, Ruppeiner and Quevedo metrics leads to existence of ensemble dependency while HPEM metric yields consistent picture. For neutral solutions, it will be shown that generalization to massive gravity leads to the presence of non-zero temperature and heat capacity for vanishing horizon radius. Such behavior is not observed for linearly charged solutions while generalization to nonlinearly one recovers this property.

Motivated by recent developments of BTZ black holes and interesting results of massive gravity, we investigate massive BTZ black holes in the presence of Maxwell and Born-Infeld (BI) electrodynamics. We study geometrical properties such as type of singularity and asymptotical behavior as well as thermodynamic structure of the solutions through canonical ensemble. We show that despite the existence of massive term, obtained solutions are asymptotically (a)dS and have a curvature singularity at the origin. Then, we regard varying cosmological constant and examine the Van der Waals like behavior of the solutions in extended phase space. In addition, we employ geometrical thermodynamic approaches and show that using Weinhold, Ruppeiner and Quevedo metrics leads to existence of ensemble dependency while HPEM metric yields consistent picture. For neutral solutions, it will be shown that generalization to massive gravity leads to the presence of non-zero temperature and heat capacity for vanishing horizon radius. Such behavior is not observed for linearly charged solutions while generalization to nonlinearly one recovers this property.

I. INTRODUCTION
Three dimensional black holes were first studied by Banados, Teitelboim and Zanelli which are known as BTZ Black holes [1][2][3]. The discovery of the BTZ black holes contributed to studies that are conducted in different aspects of physics. Among them one can point out: providing simplified machinery for studying different properties of black holes such as thermodynamical ones [4][5][6][7][8][9][10], contributing to our understanding of gravitational systems and their interactions in lower dimensions [11], existence of specific relations between BTZ black holes and effective action in string theory [12][13][14][15] and possible existence of gravitational Aharonov-Bohm effect due to the noncommutative BTZ black holes [16]. In addition, several studies were conducted in context of AdS/CFT correspondence [17][18][19], quantum aspect of three dimensional gravity, entanglement, quantum entropy [20][21][22][23][24][25][26] and holographic aspects of such solutions [27][28][29]. Furthermore, three dimensional black holes and their thermodynamical properties in the presence of different matter fields with various gravity models are addressed in literature [30][31][32][33][34][35]. Motivated by the above brief discussion, we will conduct our study in three dimensional spacetime.
On the other hand, Einstein theory of gravity introduces gravitons as massless particles with two degrees of freedom, 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 as well. The first attempt for constructing a massive theory was referred to the works of Fierz and Pauli [36]. This theory was built in a flat background and its generalization to curved background leads to the existence of a typical ghost (Boulware-Deser ghost) [37]. Presence of the ghost indicates that the theory under consideration is unstable. In order to avoid such instability, several models of massive theory are proposed. One of the ghost-free massive theories was introduced in a pioneering work by Bergshoeff, Hohm and Townsend. This theory is a three dimensional massive theory and is known as new massive gravity (NMG) [38]. The construction of NMG was done based on two important facts: first, a massless graviton has no propagating degrees of freedom in three dimensions. This leads to a symmetric tensor in three dimensions with six components with a diffeomorphism symmetry. Due to this diffeomorphism symmetry, three degrees of freedom are pure gauge and the other three ones are non-dynamical. Second, a massive graviton in three dimensions has two degrees of freedom similar to a massless graviton in four dimensions. Therefore, it is possible to build a diffeomorphism invariant theory of massive gravity in which only the massive degree of freedom is present. Black hole solutions in presence of this massive gravity theory have been obtained in literature [39][40][41][42][43][44][45]. In addition, some of their stabilities have been addressed in Refs. [46,47]. As for other three dimensional massive gravities, one can regard topological massive gravity [48], supergravity extensions [49] and critical gravity [50] (for more details, see Ref. [51]).
Another class of massive gravity was proposed by de Rham, Gabadadze and Tolley (dRGT) [52,53]. Contrary to previous three dimensional theories, dRGT theory is valid in higher dimensions as well. In this theory, the mass terms are produced by consideration of a reference metric. The stability of this theory was addressed in Refs. [54,55] and it was shown that such theory enjoys absence of the Boulware-Deser ghost [54,55]. Black hole solutions of dRGT massive gravity and their thermodynamical properties have been investigated in Ref. [56]. In addition, the stability of the Schwarzschild-de Sitter black holes in presence of dRGT massive gravity has been analyzed [57] (For detailed study regarding static spherically symmetric black holes in dRGT massive gravity, we refer the reader to Ref. [58]).
Motivated by applications of gauge/gravity duality, Vegh introduced another type of massive gravity theory [59]. In this theory, similar to dRGT theory, mass terms are built by using a reference metric except that this reference metric is a singular one. Using this theory of massive gravity, Vegh has shown that graviton may behave like a lattice and exhibits a Drude peak [59]. In Ref. [60], it was pointed out that for arbitrary singular metric, this theory of massive gravity is ghost-free and stable. Charged black hole solutions in the presence of this massive gravity have been studied in Ref. [61]. Also, the existence of Van der Waals like behavior in extended phase space has been addressed in Ref. [62]. Moreover, the generalizations of such theory to include nonlinear electrodynamics, higher derivative gravity and gravity's rainbow have been done [63][64][65][66]. In addition, thermodynamic properties, phase transitions and thermal stability of such generalized models have been investigated in Refs. [63][64][65][66]. Furthermore, it was shown that the mentioned massive gravity has sensible effects in metal/superconductor phase transition. It was pointed out that for this phase transition, the maximal tunnelling current and the coherence length of junction are affected by the mass of graviton in Vegh's massive gravity model [67]. In addition, phase transition of holographic entanglement entropy in this massive gravity is investigated as well [68]. It is also notable that the existence of massive graviton results into new phenomenology for gravitational system. Among them one can name: modification in propagation's speed of the gravitational waves [69], effects on production of gravitational wave during inflation [70,71] and existence of additional polarization for gravitational waves. The mentioned arguments motivate us to study black hole solutions in the presence of Vegh's massive gravity.
In order to have a deep insight for considering three dimensional massive gravity, one may regard an electromagnetic field as a source. Maxwell theory of electrodynamics has been employed to study electrically charged systems in various branches of physics. Although Maxwell theory is a relatively successful theory in describing different phenomena, it is not flawless. One of the main problems of this theory is existence of singularity for electromagnetic field of point-like charges at their locations. In other words, the electric field of a point-like charge diverges at the charge location which leads to an infinite self-energy. In 1934, Born and Infeld (BI) [72] introduced a nonlinear electromagnetic field theory to solve infinite self-energy problem by imposing a maximum strength of the electromagnetic field. One year later, this nonlinear theory was first employed by Hoffmann in context of Einstein gravity [73]. Although BI theory has long been discarded due to its complexity, in recent two decades, different types of black holes in the presence of BI electrodynamics have been investigated [74][75][76][77][78][79][80][81][82][83][84][85].
It is worthwhile to mention that the generalization of Maxwell field to nonlinear electrodynamics leads to enriching the new solutions and introducing a new phenomenology for the system under consideration. Among various motivations for studying BI theory, one can point out obtaining BI effective action in superstring theory and absence of singularity for D-branes [86][87][88][89][90][91], and also Eddington-inspired BI theory's applications in gravitation and cosmology [92][93][94][95][96][97][98]. In addition, the functional form of BI nonlinear theory of electrodynamics has been employed as gravitational theory as well [99,100]. For duality and other aspects of BI theory, we refer the reader to some interesting papers [101][102][103][104][105][106][107]. Due to these motivations and in order to have a better picture of charged three dimensional black holes in presence of massive gravity, we study two cases of linearly (Maxwell theory) and nonlinearly (BI theory) charged black hole solutions and investigate their thermodynamical properties.
Thermodynamical aspect of black holes has been of a great interest ever since the pioneering work of Hawking and Beckenstein [109][110][111]. Recent progresses in gauge/gravity duality has also highlighted importance of black hole thermodynamics [112][113][114][115][116][117][118][119][120][121][122][123][124][125][126]. In addition, the black hole thermodynamics plays significant role in non-perturbative aspects of quantum gravity. In other words, thermodynamical aspect of black holes may provide a basis for constructing a consistent theory of quantum gravity. Considering the mentioned motivations, we study some thermodynamical properties of three dimensional charged massive black holes. Among the thermodynamical properties, thermal stability has specific importance [127][128][129][130]. In order to black holes being physical thermodynamical systems, they should be thermally stable. In case of unstable solution, system may go under phase transition to acquire stable state. In this paper, we also study thermal stability and phase transition of black holes through canonical ensemble.
Another method of studying thermodynamical structure of black holes is through the use of geometrical method. In other words, by considering a thermodynamical quantity as a potential and its corresponding extensive parameters, one can build a phase space of the thermodynamical system. The Ricci scalar of this phase space carries information regarding thermodynamical properties. The divergencies of thermodynamical Ricci scalar may coincide to bound or phase transition points. So far, four different geometrical methods were introduced; Weinhold [131,132], Ruppeiner [133,134], Quevedo [135,136] and HPEM [137][138][139]. Several studies in context of black holes with consideration of such methods have been done [140][141][142][143][144][145]. In addition, it was pointed out that one can employ geometrical thermodynamics to obtain phase transition points in superconductors [146]. Furthermore, a comparative study in context of black holes was done in Ref. [147]. Thermodynamical concepts indicate that studying a thermodynamical system should lead to same results irrespective of employed ensemble. In other words, no ensemble dependency must exist. The existence of ensemble dependency in geometrical thermodynamics [148] could be observed through two cases: I) bound and/or phase transition points is/are not matched with divergency of thermodynamical Ricci scalar. II) extra divergency exists for the Ricci scalar which is not matched with any mentioned points. In several studies, it was pointed out that employing Weinhold, Ruppeiner and Quevedo metrics leads to existence of ensemble dependency [137][138][139]. In order to avoid such problem, HPEM metric was proposed. Motivated by the applications of geometrical thermodynamics and possible existence of ensemble dependency, in this paper, we conduct a comparative study regarding different approaches of geometrical thermodynamics for three dimensional charged black holes in massive gravity. Linearly and nonlinearly charged three dimensional black hole solutions and their thermodynamical properties were investigated in Refs. [137,149,150]. In this paper, we generalize our solutions to include massive gravity. We examine the effects of this generalization on thermodynamical aspect of charged three dimensional black holes.
The outline of the paper will be as follow. In Sec. II, we introduce action and basic field equations related to three dimensional charged black holes in the massive gravity context. Section III is devoted to obtain black hole solutions of the massive gravity in linearly charged set up. The thermodynamics and thermal stability will be investigated and some remarks regarding neutral case are given. Next, generalization to nonlinear electrodynamics is done and the effects of this generalization on thermodynamics and phase transition of the black holes are investigated. In addition, it will be shown that these black holes, despite generalization of the massive gravity and nonlinear electromagnetic fields, suffer the absence of Van der Waals like behavior in their phase spaces. Furthermore, geometrical thermodynamics is employed to study phase transition of these black holes. The paper is concluded by some closing remarks.

II. BASIC EQUATIONS
Here, we start with the 3-dimensional action of Einstein-massive gravity with an abelian U (1) gauge field where R and L(F ) are, respectively, the scalar curvature and an arbitrary Lagrangian of electrodynamics, Λ is the cosmological constant and f is a fixed symmetric tensor. Also, F = F µν F µν is the Maxwell invariant, in which F µν = ∂ µ A ν − ∂ ν A µ is the Faraday tensor and A µ is the gauge potential. In addition, c i 's are some constants and U i 's are symmetric polynomials of the eigenvalues of the d × d matrix K µ ν = √ g µα f αν which can be written as follows Taking into account the action (1) and using the variational principle, we can obtain the field equations related to the gravitation and gauge fields as where and χ µν is the massive term with the following form

III. EINSTEIN-MAXWELL SOLUTIONS IN THE CONTEXT OF MASSIVE GRAVITY:
In this section, we regard linearly charged three dimensional solutions in the context of massive gravity and investigate their geometric and thermodynamic properties.

A. Linearly charged black hole solutions
Here, we obtain static charged black holes with (a)dS asymptotes. For this purpose, we consider the metric of 3-dimensional spacetime with (− + +) signature in the following explicit form where f (r) is an arbitrary function of radial coordinate. In order to obtain exact solutions, we should make a choice for the reference metric. We consider the following ansatz metric [61] where c is a positive constant. Using the metric ansatz (6), U i 's are in the following forms [61] where indicate that since we are working in three dimensions, the only contribution of massive gravity comes from U 1 .
Due to the fact that we are going to study the linearly charged BTZ solutions, we choose the Lagrangian of Maxwell field L(F ) = −F . In addition, due to our interest in electrically charged black holes in massive gravity context, we consider a radial electric field which its related gauge potential is Using the metric (Eq. (5)) with the Maxwell field equation (Eq. (3)), one finds the following differential equation where the prime and double prime are the first and the second derivatives versus r, respectively. It is easy to find the solution of Eq. (9) as where q is an integration constant which is related to the electric charge and l is an arbitrary constant with length dimension which is coming from the fact that the logarithmic arguments should be dimensionless. It is notable that the electromagnetic field tensor is F tr = q r , which does not depend on l as a measurable physical quantity. Now, we want to obtain exact solutions for the metric function f (r). For this purpose, by employing Eq. (5) with Eq. (2), we obtain the following differential equations which are corresponding to tt(or rr) and ϕϕ components of Eq. (2), respectively. After some manipulations, one can obtain the following metric function where m 0 is an integration constant which is related to the total mass of the black hole. We should note that the obtained metric function satisfies all components of the field equation (2), simultaneously. Also, it is notable that, in the absence of massive parameter (m = 0), Eq. (13) reduces to the following linearly charged BTZ black hole Now, we are in a position to examine the geometrical structure of this solution, in brief. First, we look for essential singularity(ies). The Ricci and Kretschmann scalars are which diverge at the origin and confirm that there is a curvature singularity at r = 0. For large values of radial coordinate, r −→ ∞, the Ricci and Kretschmann scalars are, respectively, 6Λ and 12Λ 2 , in which confirm that the asymptotical behavior of the solution is (a)dS for Λ > 0 (Λ < 0). In order to study the effects of free parameters on the metric function, we can plot various diagrams (see Figs. 1 and 2). These figures show that by considering specific values for different parameters, the metric function has different behaviors. The horizon properties of the charged massive BTZ black hole may be like Reissner-Nordström black holes. In other words, this black hole may have two horizons, one extreme horizon and without horizon (naked singularity) (see Figs. 1 and 2 for more details).

B. Thermodynamics
Now, we intend to calculate the conserved and thermodynamic quantities of the solutions and examine the validity of the first law of thermodynamics.
Using the definition of Hawking temperature with relation to the surface gravity on the outer horizon r + , one can find Using the Gauss's law, the electric charge, Q, can be found by calculating the flux of the electric field at infinity, Since we are working in Einstein gravity, the entropy of the black holes can be obtained by employing the area law. It is a matter of calculation to show that entropy has the following form [109][110][111][153][154][155] By using the Hamiltonian approach and/or the counterterm method, one can find the total mass of the solutions as in which by evaluating metric function on the horizon (f (r = r + ) = 0), we obtain The electric potential, U , is defined through the gauge potential in the following form Having conserved and thermodynamic quantities at hand, we are in a position to check the validity of the first law of thermodynamics. It is easy to show that by using thermodynamic quantities such as electric charge (18), entropy (19) and mass (20), with the first law of black hole thermodynamics one can 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 (21). In other word, although massive term modifies some of thermodynamic quantities, the first law of thermodynamics is still valid.
C. Thermal stability in the canonical ensemble Here, we study thermal stability criteria and the effects of different parameters on such criteria. The stability conditions in canonical ensemble are based on the sign of the heat capacity. This change of sign could happen whether when heat capacity meets root(s) or divergency(ies). The root of heat capacity (or temperature) indicates a bound point, which separates physical solutions (positive temperature) from non-physical ones (negative temperature). Whereas the roots of denominator of heat capacity (heat capacity divergencies) represent phase transition points. The negativity of heat capacity represents unstable solutions which may undergo a phase transition and acquire stable state. In order to get a better picture and enrich the results of our study, we investigate both temperature and heat capacity, simultaneously.
The heat capacity is given by the following traditional relation Considering Eqs. (17) and (19), it is a matter of calculation to show that Obtained temperature (17) consists three terms: cosmological constant, electric charge and massive parameter. The electric charge term is negative, therefore, its contribution is related to negativity of the temperature. The contribution of the massive term depends on the sign of c 1 . This term is not coupled with any order of horizon radius, which results in it being constant. For adS spacetime, Λ−term is positive and its contribution is toward positivity of temperature, the vise versa is observed for dS spacetime. In conclusion, the dS black holes with negative c 1 has negative temperature every where, which results into absence of physical black hole solutions for this case. Otherwise, the existence of physical black holes (having positive temperature) is restricted to satisfying a condition for massive term. For small values of horizon radius, the dominant term is q−term while for large values of it, the dominant term will be Λ−term. Therefore, for dS black holes, for small and large values of the horizon radius, the temperature is negative and the positive temperature may exist where massive term is dominant (for positive c 1 ). On the contrary, for adS black holes, for small values of horizon radius temperature is negative while for sufficiently large values of horizon radius, it is positive (regardless of massive term). Now, we are in a position to study bound and phase transition points of these black holes. Solving denominator and numerator of the heat capacity with respect to horizon radius leads to following solutions for phase transition and bound points, respectively, Interestingly, only for dS black holes phase transition exists and this phase transition is independent of massive term. In other words, the existence of second order phase transition is only a function of electric charge and cosmological constant. On the contrary, the bound point is observed for both dS and adS solutions. This bound point is a function of the massive parameter, Λ and electric charge. It is worthwhile to mention that for the absence of massive gravity (m = 0), only for adS spacetime, bound point exists (no real bound point for dS black holes with m = 0). This is the contribution of the massive gravity into thermodynamical structure of the charged massive BTZ black holes. For adS black holes, only one bound point is observed, whereas, for dS black holes, two bound points could be obtained. Then again, we point out that existence of two bound points for dS black holes is due to contribution of the massive gravity. It is worthwhile to mention that in obtained heat capacity, massive term is coupled with horizon radius. To finalize the study in this section, we present various diagrams for temperature and heat capacity in case of variation of different parameters for adS and dS spacetimes (see Figs. 3 -5).
It is evident that for adS spacetime, bound point is an increasing function of the electric charge (right panel of Fig.  3) and a decreasing function of massive parameter (left panel of Fig. 3) and cosmological constant (middle panel of Fig. 3). For the absence of electric charge, no bound point is observed and solutions are always thermally stable.
For dS solutions, a second order phase transition exists (a divergency for heat capacity is observed). For previous values of different parameters, we only see non-physical solutions (see Fig. 4). By suitable choices, one can find physical solutions (see Fig. 5). It is evident that having positive temperature, hence, physical solutions is a function of the massive parameter, cosmological constant and electric charge. Increasing massive terms leads to formation of a physical region with two bound points. Smaller bound point is a decreasing function of massive parameter while the larger bound point is an increasing function of it (left panel of Fig. 5). On the contrary, increasing cosmological constant and electric charge lead to absence of physical solutions. In other words, smaller bound point is an increasing function of the cosmological constat and electric charge whereas, the larger bound point is a decreasing function of them (see middle and right panels of Fig. 5). For physical solutions, there is a phase transition of larger to smaller black holes. In other words, even in physical region, only for specific range of horizon radius thermally stable solutions exist. Once more, we emphasize that existence of stable physical solutions for dS spacetime is only due to contribution of the massive gravity.

D. Neutral Massive Black Holes
This subsection is devoted to the absence of the electric charge, since one of the most interesting effects of the massive gravity and its interesting phenomenology could be better observed in this case. Considering the case of the absence of electric charge (q = 0) and taking a look at the obtained temperature and heat capacity, one can find The mentioned quantities confirm that both T and C Q are positive for Λ < 0 (adS solutions) while in order to have physical dS solutions (Λ > 0) with positive temperature, one should set r + < m 2 cc1 2Λ . This limitation on r + leads to obtain negative C Q , and therefore, we cannot obtain physically stable dS black hole. Hereafter in this subsection, we consider adS solutions for more discussions regarding vanishing event horizon.
Considering Eq. (28), one can find that for the limit r + −→ 0, temperature of the black holes will not vanish. In other words, for vanishing the horizon radius of the black holes, the Λ−term of temperature vanishes, whereas the massive term does not vanish. Interestingly, the same behavior could be observed for heat capacity as well. On the contrary, for r + = 0, total mass and entropy of black holes will vanish.
Regarding black holes evaporation, we know that the mass and horizon radius decrease during the evaporation process. One may regard the case r + = 0 as final state of black holes evaporation. The non-zero value of the temperature for evaporation indicates that black holes leave a trace of their existences even after complete evaporation. Therefore, all the information, regarding the existence of black holes, is not actually lost in evaporation. In fact, there is a contribution to temperature of the background spacetime in place of the formation of black holes. The trace of the existence of the black holes may be detected by this contribution to temperature of the background spacetime. In other words, by measuring the fluctuation of the temperature, one can determine formation and evaporation of a black hole (a factor of m 2 cc1 4π should be observed). On the other hand, considering that thermodynamical temperature of the black holes has geometrical interpretation as well. The non-zero temperature after evaporation, indicates that information regarding the existence of the black holes is stored and encoded in geometrical structure of the spacetime. Meaning that, formation and existence of black holes altered spacetime structure in a way that information regarding its existence is preserved. This information (alteration in curvature and geometrical structure of the spacetime) presents itself as a variation of the temperature in specific place. Therefore, this variation of the temperature has no external source. Its source is in the structure of spacetime itself.

IV. EINSTEIN-BORN-INFELD SOLUTIONS IN THE CONTEXT OF MASSIVE GRAVITY:
In this section, we obtain nonlinearly charged three dimensional solutions through Born-Infeld theory in the context of massive gravity. Then, we study their geometric and thermodynamic properties.

A. Nonlinearly charged black hole solutions
At first, we obtain Einstein-BI static black hole solutions in massive gravity with (a)dS asymptotes. Using the ansatz metric (Eq. (6)), the U i 's are similar to Eq. (7). Now, we consider BI Lagrangian as matter field where β is nonlinearity parameter. Considering the previous relation of gauge potential ansatz (Eq. (8)), and using Eqs. (5), (3) and (30), one finds with the following solution where Γ = 1 + q 2 r 2 β 2 . Also, the nonzero components of electromagnetic field tensor are F tr = −F rt = q rΓ . Now, we use the introduced metric (Eq. 5) and field equation of Eq. (2) to obtain metric function f (r). So, respectively, tt (or rr) and ϕϕ components are in the following forms By employing Eqs. (33) and (34), one finds It is easy to show that the obtained metric function (35), satisfy all components of Eq. (2), simultaneously. In addition, such metric function reduces to Eq. (13) for β → ∞ (or large values of r). On the other hand, in the absence of massive parameter (m = 0), the solution (35) reduces to obtained black hole solution in [149] Here, we investigate the geometrical structure of this solution. For this purpose, we look for the essential singularity(ies). The Ricci and the Kretschmann scalars are Therefore at r −→ 0 we have which confirm that there is a curvature singularity at r = 0. Also, the Ricci and Kretschmann scalars are 6Λ and 12Λ 2 at r −→ ∞ and hence, like linear solutions, the asymptotical behavior of this solution is (a)dS for Λ > 0 (Λ < 0 ). By considering specific values for different parameters, the metric function has different behaviors. The obtained black holes may behave like Reissner-Nordström black holes, and hence, the solution may be black hole with two horizons, black hole with one extreme horizon and naked singularity (see Fig. 6 for more details). In addition, this nonlinear solution has additional properties. Following the approach of Ref. [149], one finds a critical value for the nonlinearity parameter (β c ), in which for β < β c , the nonlinear solution behaves like Schwarzschild black hole (singularity cover with a non-extreme horizon).

B. Thermodynamics
Here, we calculate the conserved and thermodynamic quantities of the nonlinearly charged black hole solution and then we check the validity of the first law of thermodynamics.
Using the definition of surface gravity on the outer horizon r + (Hawking temperature), one finds where Γ + = 1 + q 2 r 2 + β 2 . Calculations regarding the total electric charge, entropy and total mass show that the charge, the entropy and the mass of this black hole are similar to Eqs. (18)- (20), respectively. Evaluating metric function on the horizon (f (r = r + ) = 0), one can obtain total mass as a function of different parameters as In addition, the electric potential, U , is in the following form Using obtained electric charge and potential with entropy, temperature and total mass, it is straightforward to examine the validity of the first law of black hole thermodynamics as It is worthwhile to mention that although nonlinearity and massive gravity modify some of thermodynamic quantities, the first law remains valid.

C. Heat capacity and thermal stability
In this subsection, we investigate the effects of BI generalization of BTZ-massive black hole on thermodynamical stability through heat capacity. Considering Eq. (24) with entropy (19) and temperature (42) of BI black hole, one can find Solving numerator and denominator of the obtained heat capacity, one can find following relations for bound and phase transition points, respectively It is evident that bound point is a function of cosmological constant, massive and BI parameters. For adS black holes, only one bound point could be observed while for dS black holes (by satisfying specific condition) two bound points could be obtained. It should be pointed out that existence of two bound points for dS configuration is due to contribution of the massive gravity. In other words, in absence of massive gravity, only one bound point is obtainable for these black holes in dS spacetime. Contrary to bound point, the phase transition point is independent of the massive gravity. No contribution of the massive gravity is observed for r c . The dS black holes enjoy the existence of second order phase transition in their phase space, while such transition does not exist for adS black holes. The phase transition point is a function of the electric charge, BI parameter and Λ. It is evident that generalization to BI nonlinear electrodynamics has affected thermodynamical structure of the solutions (dependency of the bound point and phase transition point on BI parameter). In order to elaborate the effects of different parameters, we plot some diagrams for dS and adS spacetimes (Figs. 7, 8 for adS; Figs. 9 and 10 for dS).
Plotted diagrams of adS solutions show that bound point is a decreasing function of the massive parameter (left panel of Fig. 7) and cosmological constant (right panel of Fig. 7), while it is an increasing function of the electric charge (left panel of Fig. 8) and BI parameter (right panel of Fig. 8). For the absence of electric charge, solutions are physical and stable for every horizon radius. One of the interesting results of the plotted diagrams is regarding the variation of BI parameter. For small values of this parameter, the behavior of the temperature is Schwarzschild like [149]. Here, similar to the absence of electric charge, for sufficiently small values of nonlinearity parameter, black holes are physical and thermally stable for all the values of the horizon radius. In addition, similar to neutral case, for the limit r + −→ 0 (black holes evaporation), the temperature and heat capacity are non-zero which has similar phenomenology that was stated before with one important exception; in this case black holes are not neutral. They are nonlinearly charged.
The case of dS solutions have more variations comparing to adS case. It is evident that depending on choices of different parameters, the temperature could be only a decreasing function of the horizon radius or it may acquire a maximum (see Figs. 9 and 10 for more details). The maximum of temperature is a decreasing function of the electric charge and BI parameter (Fig. 10), while it is an increasing function of the massive gravity ( Fig. 9 left panel). In place of maximum, if it is positive, there is a phase transition of the larger to smaller black holes.
For the case of temperature being a decreasing function of the horizon radius, there may exist a bound point. If bound point exists, a region of positive temperature exists. In physical region heat capacity is negative, therefore, physical solutions are instable, vise versa is observed for the absence of bound point. It is evident that for dS black holes, existence of non-zero temperature for r + = 0 is achieved here much more easier comparing to adS solutions. But in this case for r + = 0 while temperature is positive, heat capacity is negative (Fig. 9 right panel).
The basic motivation for generalizing to nonlinearity is for solving shortcomings of the Maxwell theory and introducing new phenomenology. We see that generalization to nonlinearity leads to existence of the non-zero temperature in case of charged black holes evaporation in the presence of massive gravity. Existence of the such non-zero temperature highlights the effect of the BI electromagnetic field. To summarize, one can state that generalization to massive gravity leads to existence of non-zero temperature in case of black holes evaporation for neutral black holes and generalization to nonlinear electromagnetic field leads to existence of such property for nonlinearly charged BTZ-massive black holes.
Our final study in this section is making a simple comparison between neutral, Maxwell and BI theories (see Fig.  11). It is evident that the largest bound point belongs to Maxwell theory (dotted line in Fig. 11). No bound point is observed for neutral case (continuous line in Fig. 11) and the bound points of BI theory for different β are located between these two limits (dashed and dotted-dashed lines in Fig. 11). These modifications highlight the effects of charged and nonlinearly charged generalizations. To conclude, one can point out that the physical stable black holes are formed in smaller horizon radius in BI theory comparing to Maxwell one.

V. CRITICAL BEHAVIOR
In this section, we will investigate the possible existence of Van der Waals like behavior for obtained solutions by using the proportionality between thermodynamical pressure and negative cosmological constant (adS solution). It was shown that by considering cosmological constant as thermodynamical pressure, it is possible to enrich thermodynamical structure of black holes and introduce new phenomena such Van der Waals like phase transition. The proportionality is given by Here, instead of using usual procedure to study the possibility of the critical behavior for these black holes, we employ another method which was introduced in Ref. [151]. This method is based on using the denominator of the heat capacity. By replacing cosmological constant with its corresponding pressure and solving the denominator of the heat capacity with respect to pressure, one can obtain a relation for P . This relation is different form equation of state. The maximums of the obtained relation for pressure are places in which phase transition takes place and system has Van der Waals like behavior. In other words, the pressure and horizon radius of the maximum are critical pressure and horizon radius. This method was employed in several paper [63,64,66,152]. Now, we employ this method to study critical behavior of the system.
Using the obtained heat capacities for Maxwell and BI cases (Eqs. (25) and (46), respectively) with Eq. (49), one can find following relations by solving the denominator of these heat capacities with respect to pressure, Since electric charge and BI parameter are positive quantities, the obtained relations indicate that no positive maximum for pressure exists for Maxwell and BI cases. Therefore, even by considering the generalization of the massive gravity, Van der Waals like behavior does not exists for charged three dimensional massive black holes in presence of cosmological constant. This result is expected. By taking a closer look at the heat capacity, one can see that denominator of the heat capacity is independent of the massive gravity. This independency results into absence of the Van der Waals like behavior (previously, it was shown that BTZ black holes suffer the absence of Van der Waals like behavior in their phase diagrams). A question may rise regarding consideration of the positive cosmological constant proportional to thermodynamical pressure. Calculations for this case will lead to following equations for pressure The obtained relations for pressure confirm that there is no maximum for pressure in this case. Therefore, following previously mentioned conditions for existence of phase transition, one can conclude that no phase transition exists for this case either.

VI. GEOMETRICAL THERMODYNAMICS
Here, we regard the geometrical approach toward studying thermodynamical structure of the solutions. In this approach, thermodynamical phase space of a black hole is constructed by employing one of the thermodynamical quantities as potential. Other extensive thermodynamical quantities are components of the phase space. The thermodynamical information of the system is encoded within Ricci scalar of this thermodynamical metric. The bound points and phase transitions that were obtained in studying heat capacity, are matched with divergencies of thermodynamical Ricci scalar of the constructed phase space. In other words, divergencies of the Ricci scalar are points in which system acquire bound points or meet phase transitions. Since, this is another method for studying thermodynamical behavior of the system, the results of it must be consistent with those obtained through other methods. Meaning, that ensemble dependency must not be observed.
A successful method of geometrical thermodynamics covers all bound and phase transition points that were observed in heat capacity. In other words, no mismatch between divergencies of the Ricci scalar and mentioned points should be observed and no extra divergency should exist. There are several methods toward geometrical thermodynamics. The well known ones are: Weinhold [131,132], Ruppeiner [133,134], Quevedo [135,136] and HPEM [137][138][139]. Their corresponding metrics are in following forms where M X = ∂M/∂X and M XX = ∂ 2 M/∂X 2 . Using these metrics, it was shown that denominators of Ricci scalar of these phase spaces are [137] The factors M SS and M S in denominators of the HPEM and Quevedo metrics ensure that divergencies of the Ricci scalars are matched with bound and phase transition points, whereas, the factor M QQ in Quevedo metrics provides an extra divergency for Ricci scalar of these metrics which is not matched with any corresponding point in the heat capacity. Analytical calculations show that for Maxwell and BI cases of this paper, one can find following roots for the factor M QQ r 0−Maxwell = 4 ln r + l 2 , r 0−BI = 4 2r 2 + β 2 (1 − Γ + ) − q 2 Γ + ln (1+Γ+)πr+ 2l + (1 + Γ + ) 2r 2 + β 2 ln π − q 2 + q 2 (2 + Γ + ) ln π 2 r 4 + β 4 Γ 2 which means that for Quevedo metrics, at least, one extra divergency for their Ricci scalar is observed that is not coincidence with any bound and phase transition points. Therefore, employing Quevedo metrics results into a case of ensemble dependency. As for Ruppeiner and Weinhold, it was not possible to obtain divergencies of their Ricci scalar analytically. Therefore, we use numerical evaluation and present results in some diagrams (see Figs. [12][13][14]. It is evident that divergencies of the Ricci scalar of Weinhold (left panels of Figs 12 and 14) and Ruppeiner (right panels of Figs 12 and 14) are not matched completely with bound and phase transition points. This result indicates that using these thermodynamical metrics also leads to ensemble dependency. On the contrary, the HPEM metric provides divergencies in its Ricci scalar which are coincidence with bound and phase transition points (see Figs 13 and 15). In other words, no extra divergency and mismatched are observed. Therefore, using this method results into a uniform picture regarding bound points and phase transitions. In other words, no ensemble dependency exists for this method.

VII. CONCLUSIONS
In this paper, we have considered (non)linearly charged BTZ black holes in the massive gravity context. After investigating geometrical properties, thermodynamical structure of the solutions has been investigated through various methods. It was shown that thermodynamical structure of the solutions is highly sensitive to consideration of the adS or dS spacetime. The number of bound points and phase transitions was a function of cosmological constant and was sensitive to the sign of Λ.
Interestingly, it was shown that due to contribution of the massive gravity, number of the bound points was modified while the divergencies (phase transition) of the heat capacity was independence of the massive gravity. One of the most important results of this paper was existence of non-zero temperature for vanishing horizon radius (final state of black holes after evaporation). In other words, in case of neutral black holes, for black holes evaporation, the temperature and heat capacity were non-zero, indicating that after black holes evaporation, all information of the black holes is not lost. A trace of existence of black holes is left in terms of a fluctuation of the temperature. Such property was not observed for linearly charged black holes, whereas, interestingly, generalization to nonlinear electromagnetic field recovered this property. It means that after evaporation of the nonlinearly charged BTZ-massive black holes, the trace of their existence will remain (non-zero temperature). In addition, it was shown that generalization to nonlinearity also modified the place of bound points, hence the formation of the physical black holes.
Next, it was shown that despite the generalization to massive gravity and nonlinear electromagnetic field, these black holes suffer the absence of Van der Waals like behavior in their phase diagrams.
In addition, geometrical methods was employed to study the thermodynamical structure of the solutions. It was analytically shown that due to structure of the Quevedo metrics, employing these metrics, leads to existence of ensemble dependency. In other words, the bound points and phase transitions were not matched with divergencies of the Ricci scalar of Quevedo metrics, completely and/or extra divergencies were observed. In addition, it was pointed out that Weinhold and Ruppeiner metrics also leads to presence of ensemble dependency. Whereas, the HPEM method, was free of mentioned problems and provided a uniform picture for studying thermodynamical structure of the solutions.