A Critical Assessment of Black Hole Solutions With a Linear Term in Their Redshift Function

Different theories of gravity can admit the same black hole solution, but the parameters usually have different physical interpretations. In this work we study in depth the linear term $\beta r$ in the redshift function of black holes, which arises in conformal gravity, de Rham-Gabadadze-Tolley (dRGT) massive gravity, $f(R)$ gravity (as approximate solution) and general relativity. Geometrically we quantify the parameter $\beta$ in terms of the curvature invariants. Astrophysically we found that $\beta$ can be expressed in terms of the cosmological constant, the photon orbit radius and the innermost stable circular orbit (ISCO) radius. The metric degeneracy can be broken once black hole thermodynamics is taken into account. Notably, we show that under Hawking evaporation, different physical theories with the same black hole solution (at the level of the metric) can lead to black hole remnants with different values of their physical masses with direct consequences on their viability as dark matter candidates. In particular, the mass of the graviton in massive gravity can be expressed in terms of the cosmological constant and of the formation epoch of the remnant. Furthermore the upper bound of remnant mass can be estimated to be around $0.5 \times 10^{27}$ kg.


I. INTRODUCTION
The detection of gravitational wave signals by the LIGO collaboration [1], and the measurements of the size of the shadow of the black hole Sagittarius A* by the Event Horizon Telescope [2,3], have opened new incredible opportunities for testing gravitational theories in their strong field regime. For example restricting our attention to the latter, predictions about the shadow of a Kerr-Newman black hole with [4] or without [5] a cosmological constant, or placed inside an expanding universe [6] have been formulated. From a more fundamental point of view, the observational features of the black hole shadow may constitute a valuable tool for assessing critically the very core foundation of gravitational theories like the no-hair theorem [7,8], braneworld [9], and Gauss-Bonnet corrections [10], and may also be used as indirect evidence of plasma [11], or of cosmic fluid like dark matter in the surrounding of the black hole [12]. However, testing the Kerr hypothesis does not seem an easy task because some other parameters, like scalar fields or an electric charge, may mimic the effects of the black hole spin [13][14][15]. More seriously, different gravitational theories can be consistent with the same shadow properties because these theoretical studies rely on the solution of null geodesics about the black hole which depend only on the metric tensor and of its derivatives. In fact, the same spacetime manifold can arise in many different gravitational theories: see e.g. [16] for the Kerr black hole. The same problem arises when inferring the black hole properties from the measurements of the size of its accretion disk which is determined by solving the timelike geodesics motion [17,18]. Despite such degeneracies, parameters in the metric of a certain theory might have a different physical interpretation in another theory. As a consequence, some physical processes (e.g. thermodynamics) may, in principle, allow us to distinguish the underlying theories.
A static and spherically symmetric black hole solution with a linearly increasing term (as function of the coordinate radial distance) in the redshift function is an example of degeneracy between (at least 1 ) four different gravitational theories. In fact, the spacetime metric ds 2 = g ab dx a dx b = e a · e b = −f (r)dt 2 + dr 2 f (r) + r 2 dθ 2 + r 2 sin 2 θdφ 2 , is a solution in conformal gravity [20][21][22], de Rham-Gabadadze-Tolley (dRGT) massive gravity [23], (approximate) f (R) gravity [24,25], and general relativity (in this latter case this spacetime is known under the name of generalized Kiselev black hole) [26]. In each of these theories, the parameters M , β, and Λ, which describe the same manifold at a geometrical level, carry different physical meanings. The degeneracy of the black hole metric (1) between the conformal and the massive gravity theories has been already addressed in [27], but the comparison with the f (R) gravity and general relativity frameworks has not been explored so far. The parameter β, which makes this solution different from the Schwarzschild-(anti-)de Sitter one, leads to a modified Newtonian dynamics (MOND) [28] which has been investigated in regard of deviations in the trajectories of massive bodies in the solar system which would follow from an entropic force [29], of possible phase transitions [30], and as a solution of the galaxy rotation curves problem [21]. Studies of the quasi-normal frequencies of vibrations of this black hole [31,32], and of their anomalous behavior [33], have showed that the strong cosmic censorship conjecture can be violated for some choices of the black hole parameters [34].
In this paper, we will critically assess this degeneracy by firstly describing the geometrical meaning of the parameters M , Λ, and β in terms of some curvature invariant quantities. This will allow us to provide an algorithm for a local measurements of such parameters which is independent of the previously mentioned different physical interpretations. We will focus our attention on the latter parameter and construct an appropriate combination of curvature invariants for quantifying the "distance" between our MOND solution and the Schwarzschild-(anti-)de Sitter spacetime. Then, we will illustrate the physical importance of the parameter β in the modeling of the innermost stable circular orbit (which can be taken roughly to correspond to the location of the accretion disk), and of the photon sphere (which may be used as an approximation of the shadow size). Also in this case there is a degeneracy between the different physical theories because the analysis is based on the study of the geodesic motion which is a purely geometrical concept. We will then break this degeneracy by looking at the existence of a possible remnant of the black hole (1) after evaporation. In fact, although the radius of the remnant would have the same functional dependence on the black hole parameters, the physical mass and entropy of the remnant would differ between the various gravitational models. Physically, this means that the physical mass and entropy of the remnant would depend differently on the location of the accretion disk and shadow of the progenitor. As a consequence the formation time of the remnant would be different in different physical frameworks, whose applicability in accounting for dark matter would be directly affected; we will discuss this point in particular comparing and contrasting the claims that the black hole (1) arises in a spacetime filled with massive gravitons or in a spacetime supported by an anisotropic fluid.
Our paper is organized as follows: in Sect. II we will investigate the curvature structure of the black hole (1); this section is divided into two parts: in II A we will present some curvature syzygys both in terms of scalar polynomial curvature invariants and of Cartan curvature invariants; in II B we will propose a local procedure for a measurement of the black hole parameters via the curvature invariants for taming the teleological nature of black hole spacetimes. We will continue with the analyses of the location of the accretion disk and of the shadow in Sect. III, in which we also point out how the parameter β may play the same role as the spin of a rotating Kerr black hole. The physical degeneracies explored in these two sections will be broken in Sect. IV by showing that the different physical frameworks for which (1) is a solution predict different characteristics of the black hole remnant after evaporation. We will also deepen the analysis of the roles of massive gravitons vs. anisotropic fluid in the possibility of invoking remnants as part of the dark matter budget in IV A. Finally we conclude in Sect. V by providing some physical motivations in support of our mathematical results.

II. GEOMETRICAL INTERPRETATION OF β: CURVATURE STRUCTURE OF THE BLACK HOLE SPACETIME
In this section, we will explore the geometrical meaning of the parameter β which provides a linear modification to the black hole redshift function beyond the Schwarzschild-(anti-)de Sitter solution. First of all, we will show how this parameter affects some syzygys, that is, some algebraic constraints between the curvature invariants at any generic spacetime point, and which therefore characterize the spacetime (1) [35]. Mathematically, these relationships will provide a local and invariant description of the specific black hole solution we are interested in, and this study would be a preliminary step allowing us to isolate the black hole parameters β, M , and Λ as appropriate algebraic combinations of curvature invariants. Thus, a local procedure for the measurements of the values of the black hole parameters will be provided taming the well-known problems of the teleological nature of black hole spacetimes. We remark that the curvature objects are simply a geometrical property of the manifold which are fully determined once the metric tensor (1) is provided without the need of knowing the gravitational theory in which it was discovered.

A. Curvature syzygys
In the study of the algebraic relationships between various curvature invariants characterizing the black hole spacetime (1), we will adopt three different methods, one relying on the so-called scalar polynomial curvature invariants, a second one based on the Cartan curvature invariants, and a third one exploiting the Newman-Penrose curvature scalars. We remind the reader that the Cartan curvature invariants are the components of the curvature tensors and of their derivatives computed in the canonical frame. We will also establish some quantitative relationships between the curvature quantities delivered by these different algorithms. As a side result, we will identify two appropriate curvature invariants which can separately detect locally both the black hole and the cosmological horizons through an algebraic equation.
We begin by analyzing the curvature structure of the black hole spacetime of interest in terms of the scalar polynomial curvature invariants. Let R abcd , R ab , C abcd and G ab denote the Riemann curvature, the Ricci curvature, the Weyl curvature, and the Einstein curvature tensors in the coordinate basis, respectively; let also a semicolon denote a covariant derivative. Let us introduce the following set of scalar polynomial curvature invariants 2 [36]: First of all, we should observe that can be zero in regions other than the horizon if and only if β = 0, i.e. the MOND effects are suppressed. By direct inspection, we see that the following syzygy holds: Let us now introduce the following null coframe: where i 2 = −1, in terms of which the spacetime metric (1) reads as where an overbar denotes a complex conjugation, and round parentheses stand for symmetrization. We note that the relationships l a l a = n a n a = m a m a =m am a = 0 and −l a n a = 1 = m am a hold. In the coframe (12) we obtain the following results for the non-zero Newman-Penrose curvature scalars 3 [37]: in which we have introduced as auxiliary quantity the tracefree tensor Thus, the 1-forms (12) constitute the canonical coframe for the metric tensor (1). By this we mean that while the metric tensor is invariant under the transformations constituting the Lorentz group, i.e. spin-boost and null rotations governed by and respectively, where are arbitrary functions of the manifold coordinates, the components of the curvature tensors (and thus the Newman-Penrose scalars) may change. However, once we have fixed the null coframe (12), the curvature tensor has been reduced to its canonical form for a spacetime whose symmetry group is SO(1, 3) (see appendix B in [38]). We discover that the following two syzygys between scalar polynomial curvature invariants and Newman-Penrose scalars hold: which can be combined into: This latter quantity may be interpreted as an invariant measure of the distance between the black hole solution (1) and the Schwarzschild-(anti-)de Sitter spacetime, because it vanishes if and only if the β parameter is zero. Thus, this specific combination of the scalar polynomial curvature invariants (or equivalently of the Newman-Penrose scalars) plays the same role of the "Kerness" and "massiveness" proposed respectively in [39] and in [40] for quantifying respectively the distance between axisymmetric Kerr black holes and more realistic distorted astrophysical black holes, and analogously the distance between Banados-Teitelboim-Zanelli black holes arising in general relativity and in massive gravity. Although we have shown that for the black hole (1) such a distance could be constructed equivalently in terms of either the scalar polynomial curvature invariants or using the Newman-Penrose scalars, we will soon explain how these two different paths exhibit different computational applicabilities when trying to cast the black hole parameters as appropriate combinations of curvature quantities. In fact, the former method would require us to compute second degree quantities in the derivatives of the curvature components, while the latter procedure would be based on first degree curvature quantities, making it computationally more efficient. We remark that the Newman-Penrose scalar Φ 11 is non-zero even though the spacetime (1) may be regarded as a vacuum solution for the f (R) theory, for the Weyl gravity theory, and for the massive gravity theory: this is indeed in agreement with the claim that a modification of the gravity sector behaves effectively as an extra scalar field (i.e. an extra matter component entering the stress-energy tensor) in the general relativity picture [41]. The result obtained for Φ 11 is consistent also with the observation that the f (R) gravity model in which (1) was discovered reduces to general relativity on large distances at which Φ 11 → 0.
With respect to the coframe (12), we can compute the following interesting Cartan curvature invariants [42]: Thus, comparing (23) against (20), and J 4 − J 3 versus (21), two new syzygys relating scalar polynomial curvature invariants and Cartan invariants are discovered. As done before, we can re-formulate the measure of the distance between the black hole (1) and the Schwarzschild-(anti-)de Sitter spacetime as: and again, even though there are no conceptual differences between these various routes which deliver the same final result, we will clarify some important different applicabilities of these different methods in light of practical applications in numerical relativity. The two invariant curvature quantities (3) and (24) separately detect the horizon(s) of the black hole (1) because they vanish in their correspondence and only there (they do not provide any false positives). From the conceptual point of view, the meaning of this finding is that it is not necessary to know the entire evolution of the black hole spacetime, e.g. its behavior at null spatial infinity, for locating its horizon(s) which can therefore be observed by experimental devices of finite sizes, taming its teleological nature pointed out in [43,44]. Therefore, we have explicitly shown that the horizon constitutes a local property of the manifold (in agreement with the core principles of any relativistic field theory), and that a curvature invariant can be constructed for its detection once the geometrical symmetries of the spacetime are known regardless the gravitational theory behind it: this is consistent with the geometric horizon conjecture [40,[45][46][47][48][49][50][51][52][53]. On the other hand, our results are important also from the practical point of view in light of the so-called excision technique in numerical relativity: the black hole horizon constitutes a causal boundary separating the evolutions of phenomena occurring outside it from what it may happen inside; thus, the spacetime region delimited by the horizon must be removed (or excised) when performing numerical simulations of the evolution of a black hole. Standard excision techniques rely on the solution of the Raychaudhuri equation which accounts for the propagation of a bundle of light beams exploiting the change of their focusing properties taking place on the horizon [54][55][56]. This method requires integrating a differential equation with the need of knowing information on more than one spacetime point. However, according to our two proposals for finding this region it would be necessary only to search for the roots of an algebraic equation which imposes either the curvature quantity (3) or (24) to vanish making it computationally more advantageous. Moreover, the latter choice based on a Cartan invariant seems to be even more convenient than the former because it is based on a first degree quantity in the curvature rather than on a second degree one. Furthermore, in the latter approach the position of the horizon is found by looking for the zeroes of the frame derivative of the Weyl curvature; this can be related to searching for local extrema in the radial evolution of the tidal force generated by the black hole, as considered in a number of specific solutions [57,58], but our approach is grounded on the geometric horizon conjecture rather than on the solution of the geodesic deviation equation.

B. A local measurement of the black hole parameters
So far we have explained why some appropriate curvature quantities are computationally useful for detecting the horizon of the black hole (1). However, numerical simulations of black hole systems require to know also the values of the parameters characterizing the particular solution in hand, which are M , Λ, and β in our case. Such values are usually estimated from the size of the area of the horizon [59]; since it is necessary to evaluate an integral for computing this area, a procedure has been proposed in [39] for isolating the black hole parameters for the Kerr black hole, like its mass and spin, as appropriate algebraic combinations of some scalar polynomial curvature invariants. It is conceivable, although not proved yet, that similar results may hold for any black hole solution. As for the analysis dealing with the detection of the horizon, this finding would confirm that the physical quantities characterizing the black hole are indeed local, but it would also provide a more convenient procedure for estimating them in numerical computations. Thus, the seminal result of [39] has been extended to the case of the Kerr-Newman-NUT-(Anti)-de Sitter black hole in [51] by using both scalar polynomial curvature invariants and Cartan curvature invariants separately, and to the case of the massive Banados-Teitelboim-Zanelli black hole in [40] by using a mixture of these two types of curvature quantities.
However, these results are extremely non-trivial because a constructive procedure which can deliver the black hole parameters in terms of curvature invariants does not exist in literature. This means that the physical quantities must be isolated by hands case by case looking at the specific black hole metric under investigation through a trial-anderror procedure. Despite the difficulty of this task, we were able to relate the parameters of the black hole solution (1) to appropriate combinations of some curvature invariants by exploiting the results about its curvature structure exhibited previously. Explicitly, from the scalar polynomial curvature invariants we can get: and then the mass M can be obtained equivalently from (22) or from (28) using either scalar polynomial curvature invariants, or Newman-Penrose scalars, or Cartan curvature invariants. Moreover, using only Newman-Penrose scalars we also obtain while using only Cartan curvature invariants we can isolate the cosmological constant as While there are no conceptual differences in these alternative but equivalent algorithms for measuring locally the black hole parameters, they exhibit different applicability when practical computations are performed. In fact, we have shown that the method relying on the Newman-Penrose scalars or on the Cartan curvature invariants is preferable than the one using scalar polynomial curvature invariants, as they are less expensive computationally, being fully based only on first degree quantities. We remark also that Cartan invariants are foliation-independent [42]. It should be noted that as for the discussion about locating the horizon by using an appropriate scalar polynomial curvature invariant or a Cartan invariant, we have used here only the geometrical properties of the black hole solution (1) and not its correspondence to any specific gravitational framework.

III. PHYSICAL INTERPRETATION OF β: ACCRETION DISK AND SHADOW
In this section, we will present the physical relevance of the parameter β by claiming that its non-zero value can mimic the spin of a Kerr black hole in accounting for the size of the accretion disk of black holes. The idea that some parameters of static black hole solutions can be astrophysically indistinguishable from the spin of a rotating black hole in reference to this dataset has already been pioneered for nonsingular black holes in nonlinear electrodynamics [60], charged stringy black holes [61], deformed black holes [62], and black holes surrounded by dark matter [63]. In the case of our paper, we are in an even more extreme situation since we cannot provide an unambiguous physical interpretation to the parameter β because these analyses rely on the study of the geodesic motion which depends on the metric tensor and not the gravitational theory formulation. To our knowledge, this degeneracy has never been pointed out, despite the geodesic motion in the spacetime (1) has been extensively studied in the literature [25,[64][65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82][83].
To begin with, the Novikov-Thorne model, which is widely adopted for the description of the formation of an accretion disk around a black hole, states that the particles constituting the accretion disk follows nearly geodesics paths well approximated by the innermost stable circular orbit (ISCO) falling into the central massive object due to the effects of its gravitational field [17,18]. This makes the ISCOs especially relevant in astrophysics. In the spacetime of a static black hole, the ISCO is located at the radius r = r ISCO such that: where E = E/m is the specific energy of the test particle with mass m and energy E (which is conserved during the geodesic motion because the spacetime is static), and the effective gravitational potential should be written in terms of the black hole redshift function and of the specific angular momentum of the test particle L as V eff = f (r)(1 + L 2 /r 2 ) [84]. For the black hole (1)  , wherer := r M ,Λ := ΛM 2 ,β := βM .
Using this result and that for the ISCO in a Kerr spacetime [85], where the plus/minus signs are to be applied for test particles following retrograde/prograde orbits, we can express the degeneracy between the parameters (β, Λ) and the spin parameter a as: We remark that (34) is consistent with the asymptotically flat limiting caser ISCO = 6 whenΛ =β = 0. Furthermore, we can approximate the size of the shadow with the radius r p of the photon sphere which, for a static and spherically symmetric black hole, is given by [11,12]: Therefore, for the case in (1) we obtainr where the ± signs hold for positive/negative β respectively, and in the latter case as long the mass is such that M < −1/(6β). It should be noted that the cosmological constant does not enter this expression, and that once again this result is not sensitive to the gravitational theory adopted for interpreting the black hole metric (1). For small positive β we can approximater which shows that the photon sphere is smaller than the one in a Schwarzschild spacetime; this is because for escaping the photons need to climb a larger potential well. Physically, this may be interpreted as some amount of energy that the photons may lose when scattering against the cosmic fluid (if we adopt the Kiselev interpretation) or for winning the effects of massive gravitons. Should the size of the shadow be known from astrophysical observations, we can infer the value of the parameter β asβ Implementing this result into the information about the size of the accretion disk (34) delivers the following estimate of the rescaled cosmological constant: . (42) By inserting back the parameter M in (34), and eliminating it through the non-rescaled version of (39), we can obtain: which will be useful in what follows.

IV. BREAKING THE DEGENERACY: MASS AND ENTROPY OF THE REMNANT
The discussion about the curvature structure of the spacetime manifold, and of the properties of its accretion disk and shadow we have carried out in the previous sections are purely geometrical results. We will now show that the degeneracy between conformal gravity, massive gravity, f (R) gravity, and general relativity, which all admit the black hole solution (1) can be broken by looking at the entropy and at the physical mass of the black hole. We will also propose a scalar-tensor model which can provide as well a black hole solution with a linear term in the redshift. Our main objective is to illuminate the different physical meaning that the black hole parameters M , β, and Λ assume in each of these theoretical frameworks.
We recall that conformal gravity is derived from the action in which C αβγδ denotes the Weyl curvature tensor; in this case all the three parameters M , Λ and β do not enter the Lagrangian formulation but appear only as purely mathematical (with no physical meaning) integration constants of the field equations. Indeed in [65] it was claimed that the Modified Newtonian Dynamics (MOND) effects encoded in the linear part of the gravitational potential affecting the equations of geodesic motion do not have a well-justified physical ground because the authors interpreted the black hole solution as arising in conformal gravity. Furthermore, assuming (1) to be a solution of (44) the process of black hole evaporation was investigated in [86,87] showing that the black hole physical mass is 4 : On the other hand, in massive gravity the parameters M and Λ are still mathematical integration constants of the field equations, but instead β enters already at the Lagrangian level and it is interpreted to be (related to) the graviton mass m: this massive gravity theory is based on the action where R is the Ricci scalar and U is the modified gravitational potential. For example, in [88] the existence of a black hole remnant after evaporation was interpreted as due to the role played by a massive graviton; it was also shown that the radius, mass and the entropy of the remnant are 5 : where R r is the radius of the configuration. Two solutions for the radius of the remnant exist. It should be noted however that for Λ > 0 the solution with the minus sign should be ignored because a well-defined R r > 0 would require β > β 2 + Λ, i.e. Λ < 0. For Λ < 0, the well-posedness of R r would instead be written as β < β 2 + Λ (because a multiplication by a negative factor switches the direction of the inequality) which is fulfilled if β < − |Λ|. Moreover, again taking into account that a multiplication by a negative factor switches the sense of the inequality the condition M r > 0 reads as β 2 + Λ < 2Λ β + β. Here, the RHS is positive if β 2 < −2Λ; then we would obtain −3Λβ 2 < 4Λ 2 , i.e. 3β 2 < −4Λ should hold, which is fulfilled if −2 |Λ|/3 < β < − |Λ|. About the plus sign: for a positive cosmological constant it is guaranteed that R r > 0, while the mass of the remnant is automatically well-defined if β > 0, while should this parameter be negative, this condition becomes β 2 + Λ < − 2Λ β + β ,

Rr
Restrictions on the parameters  which is consistent because β 2 + 2Λ > 0, implying −3Λβ 2 < 4Λ 2 which holds noting that the LHS is negative. On the other hand, for a negative cosmological constant we would need to satisfy the condition β 2 + Λ < −β which can be achieved again if β < − |Λ|. Then, imposing a positive M r and repeating the same steps as previously, we obtain β 2 + Λ > − 2Λ β + β , requiring us to restrict the range β < − 2|Λ|; then we can write 3β 2 > −4Λ which is satisfied in our range. These restrictions are summarized in Table I. It should be remarked that there is no ambiguity on which sign should be considered in front of the square root in the expression for R r once the range of the parameters is decided, and that remnants in massive gravity can exist in both asymptotically de Sitter and anti-de Sitter spacetimes (the analysis in [88,89] only assumed one possibility for the sign in (48)).
Massive gravity does not belong to the class of scalar-tensor theories [90], and thus when we will show later that the black hole (1) can arise in the Brans-Dicke framework, it should be considered as a nontrivial correspondence.
Furthermore, the black hole spacetime (1) is a solution also of the (approximate) theory in which the cosmological constant Λ appears already at the level of the Lagrangian. In this formulation R c is a constant of integration, and R 0 = 6α 2 /d 2 accounts for the two free parameters α and d whose effects can be summarized into the effective black hole parameter β = α/d. The gravity model (50) has been proposed in literature because it can pass solar system tests on small length scale regime, i.e. in the limiting case r ≪ d, R ≫ Λ, and R/R 0 ≫ 2/α, in which the Lagrangian can be approximated as 6 On cosmological scales (50) reduces to the Einstein-Hilbert Lagrangian f (R) ≃ R + Λ adopted in general relativity. Therefore, the parameter β enters the Lagrangian as it does in the massive gravity theory and it quantifies again the deviations from general relativity, but is not directly related to the graviton mass. We can compute the physical mass of the black hole (1) interpreted as a solution of the theory (51), i.e. its Misner-Sharp mass which is the conserved charge associated with the timelike Killing Vector Field [92], by following the procedure of [93,94]: Here R = r denotes the areal radius, and r H the horizon location for which f (r H ) = 0. Therefore, we obtain: We remark also that any f (R) gravitational theory is equivalent to the Brans-Dicke framework where the potential V of the scalar field φ is provided by [95][96][97][98] with Φ = df (R) dR . Therefore: where we have used also that potentials which differ by an additive constant are physically equivalent. Thus, in this latter formulation β is the scaling factor of a Brans-Dicke potential. It should be emphasized at this point that the thermodynamical derivation of the remnant radius exhibited in [88] is not affected by the underlying gravitational theory. However, the remnant entropy, if now we claim that (1) is a solution of (51) would not be any longer (49) but 7 [99][100][101][102][103]: where the restrictions on the applicability of the solutions previously discussed are understood. Moreover, the manifold (1) is a solution in general relativity too, given an appropriate energy-momentum tensor. However, this should not be taken as a consequence of the correspondence between the f (R) gravity formulation in empty space and the general relativity formulation in which the spacetime is filled with an effective nonideal fluid p = ω(ρ)ρ [41]. In fact, the generalized Kiselev black hole is supported by a nonperfect fluid whose radial and tangential pressure differ from each other [104]. Therefore, the existence of a black hole remnant after the evaporation of (1) does not necessarily require the existence of a massive graviton (the thermodynamical consideration which led to the estimate of the radius of the remnant R r in [88] are unaffected if we switch to the general-relativistic interpretation), because it can be due as well to the effect of a cosmic fluid with an anisotropic pressure surrounding the black hole. The Misner-Sharp mass of (1) in general relativity is: Therefore when we set r H = R r the remnant mass would be positive whenever R r is, and thus a remnant exists in both asymptotically de Sitter and anti-de-Sitter spacetimes. We note that unlike the case of massive gravity, now we can use the plus sign in the formula for R r if Λ > 0 and if Λ < 0 along with β < − |Λ|, and the minus sign only in the latter interval. Last but not least, following [19] we will now construct a gravitational theory for the black hole (1) based on where F is an arbitrary function of its argument, by assuming F (Φ) = Φ = 1 ω(Φ) and V (Φ) = 0, and restricting to the case Λ = 0 for which an analytical result can be delivered. We need to solve [19,Eq.(6)]: under the assumptions of staticity and spherical symmetry for the scalar field Φ(r). Therefore: 7 Here we are using (9). and we obtain: where Φ 0 is an arbitrary integration constant.
A. Massive gravitons or anisotropic fluid?
The mass of the black hole remnant predicted by massive gravity (48) would be larger than the one predicted by general relativity (60) as long as which necessarily requires a positive β, i.e. to consider the plus sign in R r together with Λ > 0 (see Table I). Then, for the radius of the remnant, the previous condition can be recast as which delivers the constraint Λ < 3β 2 . Therefore, for a smaller (positive) value of the cosmological constant the mass of the remnant in massive gravity can be larger than that in general relativity; this can be achieved when the pressure of the anisotropic fluid (in the Kiselev picture) or the energy budget of the gravitons is dominating over the cosmological constant. Recalling the proposal that remnants of primordial black holes may contribute to the dark matter budget [105][106][107], our result intuitively implies that massive gravity predicts that fewer remnants are required for accounting for the same total amount of postulated dark matter. On the other hand, if we approximate the linear momentum of the remnant as p ∼ M r [108], our result would suggest also that massive gravity may provide plausible candidates for "warmer" dark matter than general relativity. Furthermore, by [109,Eq.(21)] where H form is the Hubble function at the black hole remnant formation time, we can claim that for 0 < Λ < 3β 2 the formation of the remnant occurs at a smaller H, i.e. at an earlier epoch, in massive gravity than in general relativity, while for Λ > 3β 2 the otherwise occurs. More quantitatively, if we write where γ ≃ 0.2 if the remnant has formed in the early universe [110], we will obtain in massive gravity, where we have defined Furthermore, by plugging (43) into the previous equation in the case of general relativity, we obtain the following condition on the parameter β: a 2 = 12H 2 form r 5 ISCO + 12r p γ 2 r 2 ISCO + 32H form γr 4 ISCO − (20H 2 form r 2 p + 40H form γr p + 9γ 2 )r 3 ISCO − γ 2 r 2 p r ISCO + 4γ 2 r 3 p , from which the entropy and physical mass of the remnant can be written in terms of its formation time, and of the size of the accretion disk and of the shadow of its progenitor. The entropy and physical mass estimates for the case of massive gravity can be expressed in a parametric form in terms of the same quantities by following the same procedure.
Applying the parametrization (75) to the massive gravity case in (69), and considering a series expansion for small m ≃ 0, we obtain the second-order equation from which we can then eliminate the parameter δ via the cosmological constant getting Therefore the mass of the graviton in terms of the cosmological constant and of the formation epoch of the remnant is This result is well defined if 2γ(4γ − 3H form )Λ > H 2 form , i.e. H form > 4γ/3 assuming Λ < 0 which constitutes a lower bound on the epoch of formation of the remnant. By using (68), an upper bound is set on the mass of the remnant as M r < 3/8. Restoring International System units by multiplying by c 2 /G, we get the bound on the remnant mass of M r < 0.5 · 10 27 kg, which is 1/4000 solar masses. We remark that we have never used any specific value for γ in the steps leading to this estimate.
On the other hand, if we specifically impose the linearized limit of our massive gravity paradigm to admit two degrees of freedom, the mass of the graviton should be m 2 = −2Λ/3 [111]. Then, by (80) the formation time of the remnant predicted by massive gravity is H form ≈ (−Λ + Λ 2 + 8Λ)γ (81) which is well-defined if Λ < −8. In this case by using (68) the mass of the remnant would be given by M r ≈ 1 2(−Λ+ √ Λ 2 +8Λ) .

V. CONCLUSION
In this paper, we have investigated the role that a linear term in the redshift function of a black hole is playing from both the geometrical and astrophysical perspectives. The strength of this term is quantified via the parameter β, and it arises in few different theories (general relativity, Weyl, massive and f (R) gravity). Physically the linear term has different origin. In the general relativity framework, in which our black hole spacetime is known as the generalized Kiselev solution, the black hole is surrounded by an imperfect fluid, and the linear term is due to an uneven balance between radial and transversal pressure with p r − p t = β at fixed radius. In the massive gravity interpretation of the same black hole manifold such behavior would be caused by an higher value of the graviton mass. However, observables based only on the metric and its derivatives cannot distinguish between these theories.
On the other hand, the metric degeneracy is broken when thermodynamics is taken into account. Specifically, when interpreted in light of these physically different theories the black hole entropy may be different. Furthermore the Hawking evaporation process of the black hole would deliver remnants with different values of their physical masses with direct consequences on their viability as dark matter candidates.