A class of black holes in dRGT massive gravity and their thermodynamical properties

We present an exact spherical black hole solution in de Rham, Gabadadze, and Tolley (dRGT) massive gravity for a generic choice of the parameters in the theory, and also discuss the thermodynamical and phase structure of the black hole in both the grand canonical and the canonical ensembles (for the charged case). It turns out that the dRGT black hole solution includes other known solutions to the Einstein field equations, such as the monopole-de Sitter–Schwarzschild solution with the coefficients of the third and fourth terms in the potential and the graviton mass in massive gravity naturally generates the cosmological constant and the global monopole term. Furthermore, we compute the mass, temperature and entropy of the dRGT black hole, and also perform thermodynamical stability analysis. It turns out that the presence of the graviton mass completely changes the black hole thermodynamics, and it can provide the Hawking–Page phase transition which also occurs for the charged black holes. Interestingly, the entropy of a black hole is barely affected and still obeys the standard area law. In particular, our results, in the limit mg→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_g \rightarrow 0$$\end{document}, reduced exactly to the results of general relativity.


Introduction
The question of whether a mass term for the graviton field can be introduced existed in Einstein's theory of general relativity since its inception. Massive gravity came into existence as a straightforward modification of general relativity by providing consistent interaction terms which are interpreted as a e-mail: sgghosh@gmail.com b e-mail: l_tannukij@hotmail.com c e-mail: pitbaa@gmail.com a graviton mass. Such a theory can describe our Universe, which is currently undergoing accelerating expansion without introducing a bare cosmological constant. Massive gravity modifies gravity by weakening it at a large scale compared with general relativity, which allows the Universe to accelerate, while its predictions at small scales are the same as those in general relativity. Furthermore, if a solution exists in this theory, it may also elucidate the dark energy problem. Hence, in recent years there were numerous developments in the massive gravity theories [1][2][3][4][5][6]. The first attempt was done, in 1939, by Fierz and Pauli [1]. They added the interaction terms at the linearized level of general relativity but later it was found that their theory suffered from a discontinuity in its predictions, which was pointed out by van Dam, Veltman, and Zakharov, the so-called van Dam-Veltman-Zakharov (vDVZ) discontinuity [2][3][4]. This discontinuity problem invoked further studies on the nonlinear generalization of Fierz-Pauli massive gravity. During the search for such generalizations, Vainshtein found that the origin of the vDVZ discontinuity is that the prediction made by the linearized theory cannot be trusted inside some characteristic "Vainshtein" radius and he also proposed the mechanism for the nonlinear massive gravity which can be used to recover the predictions made by general relativity [5]. At the same time, Boulware and Deser found that such nonlinear generalizations usually generate an equation of motion which has a higher derivative term yielding a ghost instability in the theory, later called a Boulware-Deser (BD) ghost [6]. However, these problems, arising in the construction of the massive gravity, have been resolved in the last decade by first introducing Stückelberg fields [7]. This permits a class of potential energies depending on the gravitational metric and an internal Minkowski metric. Furthermore, to avoid the reappear-ance of the ghost in massive gravity, the set of allowed mass terms was confined and furnished perturbatively by de Rham, Gabadadze, and Tolley (dRGT) [8,9]. They summed these terms and found three possibilities, viz. quadratic, cubic, and quartic combinations of the mass terms. The dRGT massive gravity is constructed suitably so that the equations of motion have no higher derivative term, so that the ghost field is absent. However, these nonlinear terms lead to complexity in the calculations and hence, in general, finding exact solutions in this theory is strenuous. Nevertheless, recently, several interesting measures have been taken to obtain the spherically symmetric black holes in various massive gravities [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27]. In particular, a spherically symmetric black hole solution with a Ricci flat horizon in four-dimensional massive gravity with a negative cosmological constant was obtained by Vegh [10] and was generalized to study the corresponding thermodynamical properties and phase transition structure [11][12][13]. The spherically symmetric solutions for dRGT were also addressed in [14,15], the corresponding charged black hole solution was found in [16], and its bi-gravity extension was found in [17], which includes as particular cases the previously known spherically symmetric black hole solutions.
(See [18][19][20], for reviews on black holes in massive gravity, see also [21] for the black hole solution in other classes of massive gravity.) The main purpose of this paper is to present a new class of exact spherically symmetric black hole solutions including a generalization to the charged case in dRGT massive gravity, and also to discuss their thermodynamical properties. It turns out that the solution discussed in this paper represents a generalization of the Schwarzschild solution that includes most of the known dRGT black hole solutions. The paper is structured as follows. In the next section, we review dRGT massive gravity. We present the modified equations of motion for dRGT massive gravity and a class of exact black hole solutions in Sect. 3. The calculation of the thermodynamical quantities associated with the dRGT black hole solution and the study of the phase structure of a black hole in the canonical ensemble approach are the main subject of Sect. 4. The analyses in Sect. 4 are extended for the charged case in Sect. 5, and finally we summarize our results and evoke some perspectives to end the paper in Sect. 6. We have used units which fix the speed of light and the gravitational constant via 8π G = c 4 = 1.

dRGT massive gravity
We begin by reviewing dRGT massive gravity, which is a well-known nonlinear generalization of a massive gravity and is free of the BD ghost by incorporating higher order interaction terms into the Lagrangian. dRGT massive gravity can be represented as Einstein gravity interacting with a scalar field, and hence its action is the well-known Einstein-Hilbert action plus suitable nonlinear interaction terms as given by [9] where R is the Ricci scalar and U is a potential for the graviton which modifies the gravitational sector with the parameter m g interpreted as graviton mass. Moreover, the action is written in units such that the Newtonian gravitational constant is unity (thus, κ 2 ≡ 8π ). The effective potential U in four-dimensional spacetime is given by in which α 3 and α 4 are dimensionless free parameters of the theory. The dependencies of the terms U 2 , U 3 , and U 4 on the metric g and scalar fields φ a are defined as where where f ab is a reference (or fiducial) metric and the rectangular brackets denote the traces, namely [K] = K μ μ and [K n ] = (K n ) μ μ . The four scalar fields φ a are the Stückelberg scalars, which are introduced to restore general covariance of the theory. Note that, generally, there are additional mass terms for the theory in higher-dimensional spacetime which are provided explicitly in Appendix A. One may recognize the interaction terms as symmetric polynomials of K; for a particular order, each of the coefficients of possible combinations are chosen so that these terms will not excite higher derivative terms in the equations of motion. Actually, this definition of K is not unique since it is possible to have the same action with a different definition of K-the alternating action is given in Appendix A.
To proceed, we choose the unitary gauge φ a = x μ δ a μ [10]. In this gauge, the tensor g μν is the observable metric describing the five degrees of freedom of the massive graviton. Note that since the Stückelberg scalars transform according to the coordinate transformation, once the scalars are fixed, for example, due to choosing the unitary gauge, applying a coordinate transformation will break the gauge condition and then introduce additional changes in the Stückelberg scalars. Also, we redefine the two parameters α 3 and α 4 of the graviton potential in Eq. (2) by introducing two new parameters α and β, as follows: By varying the action with respect to metric g μν , we obtain the modified Einstein field equations as where X μν is the effective energy-momentum tensor obtained by varying the potential term with respect to g μν , In addition to the modified Einstein equations, one can obtain a constraint by using the Bianchi identities as follows: where ∇ μ denotes the covariant derivative which is compatible with g μν . Henceforth, we shall use α and β, instead of the parameters α 3 and α 4 .

Black hole solution in dRGT massive gravity
In this section, we will look for a static and spherically symmetric black hole solution of the modified Einstein equations (8) with the physical metric ansatz The solution is found and classified into two branches: d(r ) = 0 or h(r ) = h 0 r where h 0 is a constant in terms of the parameters α and β [28][29][30]. The most interesting branch is the diagonal branch, d(r ) = 0, since it is simpler to analyze. For example, a class of a charged black holes within the diagonal branch in dRGT massive gravity was also investigated in Refs. [31]. The exact solutions for this ansatz are complicated and thus it is difficult to use these solutions to analyze the properties of black hole. Note that one may simplify the solution by choosing some specific relations of the parameters α and β for example α = −3β [16]. Furthermore, a class of parameters satisfying β = α 2 3 simply yields the Schwarzschild-de Sitter solution [32] (see also the analyses of such solution in Refs. [33][34][35][36]).
It is important to note that most solutions are asymptotically de Sitter or anti-de Sitter. This is not surprising since at large scale the theory should recover the cosmological solution in which the graviton mass will play the role of cosmological constant to drive the late-time acceleration of the Universe. However, a class of black hole solutions in dRGT massive gravity (or in other classes of massive gravity, e.g. the model in Ref. [37]) may encounter the issues of superluminality, the Cauchy problem, and strong coupling. (See [38] for the issues in dRGT and also [21,[39][40][41] for those in another model of massive gravity.) Since the fiducial metric seems to play the role of a Lagrange multiplier to eliminate the BD ghost, one can choose an appropriate form to simplify the calculation. In the present work, we will follow [10][11][12][13] by choosing the fiducial metric to be where c is a constant. With the choice of the fiducial metric, the action remains finite since it only contains non-negative powers of f μν (see [10], for more details). It is important to note that the effective energy-momentum tensor in Eq. (9) is derived by requiring that the fiducial metric must be non-degenerate. From Eq. (12), it is obvious that the fiducial metric is degenerate and then one may not use the effective energy-momentum tensor expressed in Eq. (9). However, as pointed out in [42], one can use the Moore-Penrose pseudoinverse of the metric K μ ν in order to find the effective energy-momentum tensor and it provides the same expression as in Eq. (9). Therefore, one can use the effective energymomentum tensor defined in Eq. (9) for the form of the fiducial metric. Moreover, the results can be checked by using the mini-superspace action as usually done in cosmological analysis. One can obtain the equation of motion by using the Euler-Lagrangian equation for the variables n and f . As a result, we found that the equations of motion still valid.
For the physical metric, we will consider the diagonal branch of the physical metric by setting d(r ) = 0. In order to obtain a black hole solution, we will choose the function h(r ) = r . Then the physical metric can be written as Up to this point, we have chosen to investigate just a class of the solutions to the dRGT massive gravity which possesses symmetries of our interest. Symmetries of solutions are of great importance in massive gravity since they affect the number of degrees of freedom and the stability of the theory at times, which are also significant issues for many massive gravity models. As found in a cosmological context, the number of degrees of freedom crucially depends on the isotropy and homogeneity of the background physical metric as well as the form of the fiducial metric [43][44][45]. Therefore, it is worthwhile to find, for this choice of the physical and the fiducial metric, whether the number of degrees of freedom still is valid and each of them is healthy. This issue is out of scope of this work and we leave this investigation for further work.
From the ansatz in Eq. (13), the components of the Einstein tensor can be written as Computing the effective energy-momentum tensor in Eq. (9) with this ansatz, the tensor X μν can be written as Note that X t t = X r r . There are specific values of the parameter c by which this effective energy-momentum tensor behaves like a cosmological constant. In particular, c = 0 simplifies each of the diagonal component of the effective energymomentum tensor so that they depend only on the parameters of the theory, which are all equal constants. Actually, by setting c = 0 in Eq. (6), the tensor K μ ν will be equal to the identity matrix leading to the fact that the mass terms are all constants at the Lagrangian level, which is corresponding to the cosmological-constant term. This is a crucially different point between our model and the one in Ref. [11]. In our model, the solution can be reduced to a Schwarzschild-AdS/dS solution where the cosmological constant-like term can be expressed in terms of the graviton mass, while the cosmological constant-like term in Ref. [11] is introduced by hand and is not related to the graviton mass.
Substituting all components of Einstein tensor and the effective energy-momentum tensor into Eq. (8), the modified Einstein equations can be written explicitly as From Eq. (20), one can obtain the solution of f which can be written as where and M is an integration constant related to the mass of the black hole. This solution incorporates the cosmologicalconstant term, namely Λ, naturally in terms of the graviton mass m g , which should not be surprising since the graviton mass serves as the cosmological constant in the selfexpanding cosmological solution in massive gravity. Moreover, this solution can be identified with the known solutions in general relativity. Note that the solutions in alternative forms with other sets of parameters are shown explicitly in Appendix A. In the case m g = 0 we have the Schwarzschild solution, as expected. For c = 0, which sets γ = ζ = 0, the solution can be classified according to the values of α and β. If (1 + α + β) < 0, the solution is in the form of Schwarzschild-de Sitter, while the case (1 + α + β) > 0 yields the Schwarzschild-Anti-de Sitter solution. Using Eqs. (20) and (21), one finds that This equation implies that the functions f and n differ only by a constant. One can choose the constant to obtain a black hole solution such that It may be verified by direct substitution of this solution into Eq. (22). The dRGT solution outlined here contains, for instance, the Schwarzschild solution (m g = 0), the de Sitter/Anti-de Sitter solutions, the global monopole solution of general relativity, and thus it also contains the monopole de Sitter-Schwarzschild solution. Note that the last term, the constant potential ζ , corresponds to the global monopole term. A global monopole solution was introduced by Barriola and Vilenkin [46], and usually comes from a topological defect in high energy physics at the early Universe resulting from gauge-symmetry breaking [47,48]. However, in this solution, the global monopole is contributed by a graviton mass. In addition, we see that the solution here is similar to the four-dimensional solution found in Ref. [11]. In Ref. [11], the author obtained the black hole solution where the bare cosmological constant is present, while, on the other hand, our result incorporates a cosmological-constant-like behavior, namely Λr 2 , in the black hole solution automatically due to introducing the potential term. However, this solution highly depends on the choice of the fiducial metric; changing to other forms of the fiducial metric will significantly affect the solution. This kind of dependency is one of the important properties of massive gravity. For example, from a cosmological point of view, one cannot have a nontrivial flat cosmological solution with a Minkowski fiducial metric [49]; only the open FLRW solution is allowed [43], where the FLRW solution with arbitrary geometry exists when the FLRW fiducial metric is considered [44]. By generalizing the form of the fiducial metric, the nontrivial cosmological solutions can be obtained [45].

Thermodynamics of the black hole
For black holes in de Sitter space, there exists more than one horizon, and the multiple horizons correspond to different thermodynamic systems. Next, we shall explore the thermodynamics of the dRGT massive gravity black hole solution given by Eq. (23) by assuming that the black hole is a closed system, i.e., there is no particle (or charge) transfer or cre-ation/annihilation. In other words, the black hole is assumed to be a canonical ensemble system. The horizons, if they exist, are given by zeros of g rr = 0 [50] or which may admit three real roots. For some classes of parameter setups, there may exist up to three horizons as shown in Fig. 1. From the left panel in this figure, we use a plot to find the region in (α, β)-space for the existence of three positive real roots of Eq. (27) by setting m g = 1, c = 1. We also pick up two points in this region to show the profile of f (r ) in the right panel of this figure. We note that the gravitational mass of a black hole is determined by f (r + ) = 0, which in terms of the outer horizon radius r + reads wherer + = r + /c. For simplicity, one can work in dimensionless parameters with units of m 2 g c 2 = 1. In this unit, the conditions for a positive value of the black hole mass can be written as [51] and hence the temperature T of the black hole becomes Taking the limit m g = 0, one recovers the temperature for the Schwarzschild black hole [51]: The crucially different point in the temperature profile compared with the one obtained in the Schwarzschild solution is that it is possible to find the positive local minimum of the temperature as shown in the right panel of Fig. 2.
In the left panel of this figure, we use units of m 2 g c 2 = 1 and we show the region matching the requirement of positivity of the local minimum of the temperature, T +(min) > 0, together with the positivity of the horizon size, r + > 0. We also adopt a simple choice of the parameters in this region, such as (α = 1, β = 0.2), to illustrate the existence of the positive local minimum of the temperature as shown in the right panel of Fig. 2. Note that, in the unit of m 2 g c 2 = 1, the dimensionless version of the temperature and horizon size can be written asT = T c andr + = r + /c. These are the actual values of the quantities plotted in the right panel of Fig. 2. For simplicity, we set c = 1, which leads to m g = 1 andT = T as well asr + = r + .
Next, we turn our attention to the important thermodynamic quantity associated with the black hole horizon which is its entropy (S). The black hole is supposed to obey the first law of thermodynamics, or dM = T dS. To calculate the entropy, we use Integrating the above equation leads to This can be written as S = A + 4 with the area of the horizon A + = 4πr 2 + , which is the famous area law. Interestingly, the graviton mass does not significantly affect the form of the entropy; it contributes only as a correction for the horizon radius which can be seen explicitly from Eq. (29). The behavior of the temperature and its corresponding horizon is shown in Fig. 2, in which there exists a local minimum of the temperature. This feature suggests that there should be a transition between two states; from the non-black hole or hot flat space state to a black hole. The transition was realized by Hawking and Page [52] (see also [53]), who found the characteristics of the so-called Hawking-Page phase transition. The transition exists if it is thermodynamically spontaneous, or alternatively, globally thermodynamically stable, as can be found by evaluating the free energies between two states. In other words, this corresponds to evaluating the Euclidean actions of those two states where the temperature is treated as a period of the imaginary time. Since in this case, there is no particle transfer, we compute the Helmholtz free energy The existence of a globally thermodynamical stability of the black hole is determined by the condition F ≤ 0. Therefore, the transition exists if where the transition from the vacuum state to the black hole takes place at F = 0. In addition to the globally thermodynamical stability, one can determine the locally thermodynamical stability by examining the sign of the heat capacity. The heat capacity of the black hole is given by By using the relations in Eqs. (28) and (32), the heat capacity becomes One can see from Fig. 3 that there is a particular horizon for which the corresponding heat capacity diverges. The divergence is due to the minimum temperature and the corresponding horizon is exactly the horizon of minimum temperature. Furthermore, requiring the locally thermodynamical stability of the black hole, the black hole must obey the condition One can see from Eqs. (37) and (40) that both of the thermodynamical stabilities depend on the parameters Λ and ζ , which are determined by the parameters of the massive gravity theory, namely, α and β. We use the plot to find the allowed region in (α, β)-space by using the unit of m 2 g c 2 = 1 and setting M = 1. We show the validity of those parameters where both stabilities are assumed in the left panel of Fig. 4. Note that one can express r + in terms of α, β, and M from Eq. (29) and then substitute this expression into the stability conditions to find stability regions in (α, β) space where M is held fixed to be a positive constant. From this figure, one can see that there exists an allowed region for the transition phase with the parameters α, β ∼ O(1). This suggests that massive gravity can naturally provide the Hawking-Page phase transition without requiring fine-tuning of the parameters. For the right panel of Fig. 4, we combine three important regions according to the previous plots including the regions satisfying thermodynamical stability conditions (intersection region in the left panel of Fig. 4), T +(min) > 0 and r + > 0 (intersection region in left panel of Fig. 2), and the existence of three horizons (region in the left panel of Fig. 1). From this figure, one can see that the thermodynamically stable region of the black hole together with positive temperature and horizon size is not compatible with the region indicating the existence of three horizons. This is due to the fact that the temperature is proportional to f (r + ) and the condition of existence of three horizons is that there exist two extremum points such that f (r ) = 0. This implies that there always exists a range of horizons (r + ) at which the black hole temperature is negative if there exist three real black hole horizons. Therefore, we use a set of parameters which give rise to one or two horizons to illustrate the thermodynamical quantities such as the temperature and the heat capacity.
It is worthwhile to note that in the expression of dimensionless variables, such as Eq. (29), the dimensionless variable of r + isr + = r + /c. Therefore, the dimensionless horizon size is inversely proportional to the parameter c. This means that the greater the value of the horizon size, the smaller the value of the parameter c. In the left panel of Fig.  5, we explicitly show the stability region in (α, β)-space with different values of the parameter c such that c = 0.8, c = 1.0, and c = 1.2. Furthermore, we also show the allowed region in (α, c)-space by fixing β = 0.1 in the right panel. From this figure, it is found that the greater value of the parameter c corresponds to the smaller value of the horizon size and the smaller allowed region in the parameter space. Therefore, the phase transition tends to occur more easily at large horizon size. Note that, in this case, the black hole is treated as a canonical ensemble system where particle transfer is prohibited. We will discuss the charged black hole case in the next section.

Charged black hole
In this section, the black hole with non-zero charge and both the grand canonical aspect of the black hole, where the charge transfer is allowed, and the canonical ensemble point of view will be discussed. Hence, it will be interesting to consider the charged generalization of the above solution. The action in Eq. (1) with Maxwell term reads where F μν ≡ ∇ μ A ν − ∇ ν A μ is the Maxwell strength tensor and A μ is a vector potential (the action is written in Gaussian units). We limit our study by considering a spher-ically symmetric dRGT black hole filled only with a static charge which is accompanied by A μ = (A(r ), 0, 0, 0). The presence of the static charge gives rise to a non-vanishing energy-momentum tensor, thus: Thus, one can rewrite the Einstein equations in Eqs. (20) and (21) in the presence of a static charge as where X t t and X r r are given in Eqs. (17) and (18). Taking these equations into account, one can find that the constraint in Eq.
where the constraint in Eq. (25) is applied in the calculation. This equation simply implies the solution of the form which is exactly an electrostatic potential in a generic electrodynamics, where k is an integration constant. The corresponding electric field is To determine the value of k, one may consider an electric field from a charge Q at large r where the spacetime is asymptotically flat which, in Gaussian units, must take the form Obviously, this implies that the integration constant k must be identified to the charge Q. By solving Eqs. (43) and (44), one finds a charged dRGT black hole solution as where Λ, γ, ζ are defined similarly to Eq. (24); note that we still have the relation in Eq. (26). In general, this solution may have up to four horizons, as illustrated in Fig. 6. Again, the innermost two horizons are the horizons that can be found in the Reissner-Nordström black hole in the general relativity theory, where the others are the cosmological horizons associated with the existence of the graviton mass. Following the same procedures as in Sect. 4, one can find the mass and the temperature evaluated at the black hole horizon as Since there is non-zero charge involved, the thermodynamics will be different. In detail, one can treat the black hole as an open system, i.e. a grand canonical ensemble, where the charge transfer is allowed while another can view the system as a closed one or a canonical ensemble where the charge is a non-zero constant.

Grand canonical ensemble
To treat it as a thermodynamical object, one can consider the black hole as a grand canonical ensemble system where the chemical potential is held fixed as μ = Q r + (in the Gaussian unit). The corresponding temperature and entropy are From Eq. (53), the temperature will be lower due to the effect of the chemical potential, μ. In the grand canonical ensemble, the corresponding free energy is given as the Gibbs free energy, From this expression, one can see that the free energy shifts to a smaller value due to the contribution from the chemical potential in the last term. In other words, the contribution from the charge makes the free energy more negative. The sign of the free energy depends on the value of 1 − Λ 3 r 2 + + ζ − μ 2 ; in other words, there exists globally thermodynamical stability when is satisfied. Moreover, one can check locally thermodynamical stability by considering the heat capacity, Requiring locally thermodynamical stability, one should have For this grand canonical ensemble aspect, both of the stabilities depend not only on the parameters Λ and ζ like in the former analysis but also on μ where the former two are again determined by α and β. In the left panel of Fig. 7, it shows the validity of those parameters where both stabilities are assumed. From this figure, one may see that it is not difficult to find the allowed region with parameters α, β ∼ O (1), which suggests that the theory naturally provides the Hawking-Page phase transition. In the right panel  Fig. 7, a comparison of the allowed regions between the charged and non-charged (μ = 0) cases is illustrated. From this plot, it is found that the existence of the charge makes the allowed region smaller. Furthermore, we can see the correspondence between the minimal temperature and the divergence of the heat capacity graphically in Fig. 8.

Canonical ensemble
On the other hand, if charge transfer is prohibited, one can consider the black hole as a closed system, or a canonical ensemble with fixed non-zero charge Q. The mass and temperature are given readily by Eqs. (51) and (52), respectively. Furthermore, the corresponding entropy still obeys the area law, Since the charge is fixed, the appropriate free energy in consideration is the Helmholtz free energy, which is Similarly, the condition for globally thermodynamical stability is Furthermore, to examine the locally thermodynamical stability, the corresponding heat capacity is computed, The condition for a locally stable black hole to exist is Here, the conditions for both stabilities depend on Λ and ζ (or, similarly, α and β), as well as on the charge Q rather than the potential μ as in the grand canonical ensemble case, as shown in the left panel of Fig. 9. Again, from the allowed region in this figure, the theory can provide the Hawking-Page phase transition naturally. From the right panel of this figure, it is also found that the existence of the charge makes the allowed region smaller, similar to the grand canonical case. This is implied by Eq. (52) where the region for which T > 0 is reduced by the presence of charge. Similarly, the correspondence between the minimal temperature and the divergence of the heat capacity can be seen in Fig. 10. The effect of charge on the stability regions of both the grand canonical ensemble and the canonical ensemble compared with the non-charged case are shown in Fig. 11. From this figure, one can see that the effect of the chemical potential μ in the grand canonical ensemble and the charge Q in the canonical ensemble decreases the size of the allowed given by the red thick curves compared with the black dashed curves corresponding to those of the neutral black hole. The divergence at large r + in the heat capacity corresponds to the minimal temperature of the black hole. Moreover, the heat capacity diverges at two values of r + region of the parameters. This can be seen from Eqs. (52) and (53), since the contribution from the charge and the chemical potential makes the positive temperature region smaller. In the canonical ensemble, we also found that the parameter region that satisfies the conditions for the existence of four horizons is not compatible with the stability region, similar to the non-charged case. The effect of the parameter c on the region of stability is also similar to the non-charged case. The stability region will increase where the horizon size increases or the parameter c decreases.

Concluding remarks
The dRGT massive gravity is a natural extension of Einstein's theory of general relativity, providing mass to the graviton, and it is a great arena for theoretical physics research. The dRGT massive gravity describes nonlinear interaction terms as a correction of the Einstein-Hilbert action and hence admits general relativity as a particular case. It is believed that dRGT massive gravity may provide a possible explanation for the accelerated expansion of the Universe that does not require any dark energy or cosmological constant. Hence, dRGT massive gravity has received significant attention [54] including searches for black holes [18]. In this paper, we have obtained a class of black hole solutions in dRGT massive gravity, and we studied the thermodynamics and phase structure of the black hole solutions. In the dRGT massive gravity outlined in this paper, there are three terms in the effective potential associated with the graviton mass. Furthermore, we note that due to the inclusion of the massive gravity term in the action, the Schwarzschild solution in general relativity is modified. Interestingly, it turns out that solutions to the Einstein field equations, such as the monopolede Sitter-Schwarzschild model, become solutions in dRGT massive gravity for suitable choices of the parameters of the theory, where the coefficients for the third and fourth terms in the potential and the graviton mass in massive gravity naturally generate the cosmological constant and the global monopole term. The corresponding thermodynamical quantities are also changed. However, the black hole entropy is not affected significantly by the existence of the graviton mass, and it still obeys the standard area law as in general relativity. Moreover, we consider the charged black hole solution in both the grand canonical and the canonical ensembles to analyze the thermodynamics and phase transition. We have demonstrated through the calculation of the heat capacity and the free energy that there is a critical point where the heat capacity diverges and a phase transition is possible without requiring fine-tuning of the parameters as shown in Fig. 4. Even though it is possible to obtain three horizons in some region of parameter space, the region is still not compatible with the stability region. This implies that the phase transition will not occur when the black hole has three horizons. The presence of the charge will not change this argument, the phase transition will not occur when the black hole has four horizons. The presence of the charge affects the appearance of the Hawking-Page phase transition such that the allowed parameter region decreases as shown in Fig. 11. We also found that the phase transition tends to occur in the large horizon size in both the charged and the non-charged cases. The black hole solutions obtained are immensely simplified due to the choice of the fiducial metric as f μν = diag(0, 0, c 2 , c 2 sin 2 θ) and the choice of the Stückelberg scalars. It will be interesting to apply the technique discussed here in other massive gravities to get black holes. It will also be interesting to consider the motion of particles in the background of the dRGT massive black holes considered and to see how the graviton mass affects the equations of motion. These and related areas are for future investigation.
where the building block tensor is Note that here the square brackets acting on the indices denote anti-symmetrization, For the four-dimensional ghost-free massive gravity, the anti-symmetric contractions of the terms that are higher than fourth order vanish, which leaves the non-zero interaction term in Eq. (A.2) as follows: To transform these terms to the different convention introduced above, one can transform the coefficients α i to c i via the transformation matrix,