Black hole thermodynamics in Lovelock gravity's rainbow with (A)dS asymptote

In this paper, we combine Lovelock gravity with gravity's rainbow to construct Lovelock gravity's rainbow. Considering the Lovelock gravity's rainbow coupled to linear and also nonlinear electromagnetic gauge fields, we present two new classes of topological black hole solutions. We compute conserved and thermodynamic quantities of these black holes (such as temperature, entropy, electric potential, charge and mass) and show that these quantities satisfy the first law of thermodynamics. In order to study the thermal stability in canonical ensemble, we calculate the heat capacity and determinant of the Hessian matrix and show in what regions there are thermally stable phases for black holes. Also, we discuss the dependence of thermodynamic behavior and thermal stability of black holes on rainbow functions. Finally, we investigate the critical behavior of black holes in the extended phase space and study their interesting properties.


I. INTRODUCTION
It is important to understand the UV behavior of general relativity, and various attempts have been made to obtain the UV completion of general relativity such that it reduces to the general relativity in the IR limit. This has been the main motivation behind the development of the Horava-Lifshitz gravity [1,2]. In the Horava-Lifshitz gravity, space and time are made to have different Lifshitz scaling. Thus, the Horava-Lifshitz gravity reduces to general relativity in the IR limit. However, its behavior in the UV limit is different from that of general relativity. Motivated by the development of the Horava-Lifshitz gravity, the UV completion of various geometric structures in the string theory has been studied by taking a different Lifshitz scaling for space and time. In fact, such a UV completion has been studied for geometries that occur in the type IIA string theory [3] and type IIB string theory [4]. The AdS/CFT correspondence has also been used for analyzing the geometries where the space and time have different Lifshitz scaling [5][6][7][8]. This formalism has also been applied for analyzing the UV completion of dilaton black branes [9,10] and dilaton black holes [11,12]. The Horava-Lifshitz gravity is based on a deformation of the usual energy-momentum dispersion relation in the UV limit, such that it reduces to the usual energy-momentum dispersion relation in the IR limit. Such a deformation of usual energy-momentum dispersion in the UV limit has been observed to occur in ghost condensation [13] and non-commutative geometry [14,15].
There is another approach for obtaining the UV completion of general relativity based on the deformation of the usual energy-momentum dispersion relation in the UV limit, and this approach is called gravity's rainbow [16]. In gravity's rainbow, the geometry of spacetime is made energy dependent and this energy dependence of the spacetime metric is incorporated through the introduction of rainbow functions. It has been suggested that the deflection of light and gravitational red-shift can be used to test various choices of rainbow functions [17].
It has been demonstrated that for a certain choice of rainbow functions, the gravity's rainbow is related to the Horava-Lifshitz gravity [18]. String theory can also be used as a motivation for gravity's rainbow. This is because the background fluxes in string theory produce a noncommutative deformation of the geometry [19,20], and noncommutativity has also been used to motivate one of the most important rainbow functions in gravity's rainbow [21,22]. In string field theory, a tachyon field can have the wrong sign for its mass squared, and the perturbative string vacuum become unstable in the presence of such a tachyon field [23]. This existence of such an unstable perturbative string vacuum spontaneously breaks the Lorentz symmetry. It may be noted that a gravitational Higgs mechanism in supergravity theories also spontaneously breaks the Lorentz symmetry [24]. The spontaneous breaking of the Lorentz symmetry deforms the usual energy-momentum relations, and this in turn can be used as a motivation for introducing gravity's rainbow.
In gravity's rainbow a one-parameter family of energy dependent orthonormal frame fields introduced. This gives rise to a one-parameter family of energy dependent metrics as, g µν (ε) = e µ a (ε)e aν (ε) in which ε = E/E p is dimensionless energy ratio, E is the energy of the test particle and E P is the Planck energy [16]. Here, the rainbow functions are used to relate these new tetrad fields to the usual frame fieldsẽ µ a of general relativity, f (ε)e µ 0 (ε) =ẽ µ 0 and g(ε)e µ i (ε) =ẽ µ i , where i is the spatial index. In the IR limit, ε → 0, we have lim ε→0 f (ε) = lim ε→0 g(ε) = 1, and so in the IR limit one recovers the general relativity. It may be noted that if the energy E was just a non-dynamical parameter in the theory then we could gauge it away by rescaling. However, it dynamically depends on the coordinates, and it originally breaks the diffeomorphisms symmetry of the full metric. In fact, even the local symmetry in gravity's rainbow is not Lorentz symmetry, as it is based on deformation of the usual energy-momentum dispersion relation, which breaks the Lorentz symmetry in the UV limit (similar to Horava-Lifshitz gravity) [16]. As the energy ε is an implicit function of the coordinates, and the the rainbow functions are explicit functions of the energy ratio ε, they are also dynamical functions of the coordinates, and so they cannot be gauged away. It may be noted that for certain systems, an explicit dependence of the energy on the coordinates has been constructed [18]. Even though it is difficult to find such an explicit dependence of rainbow functions on energy for different systems, it is important to know that these rainbow functions are implicitly dynamical functions of the coordinates, and so they cannot be gauged away. Thus, these functions are expected to produce physically different results from general relativity [25]. In this paper, we will explicitly demonstrate that the thermodynamics of the rainbow deformed Lovelock black hole is different from a black hole in the usual Lovelock theory.
It may be noted that the rainbow functions physically change the thermodynamics of black holes [25]. In fact, the rainbow deformation of a black hole produces non-trivial effects like the existence of a black hole remnant at the last stage of the evaporation of a black hole, and this has been observed to have phenomenological consequences for the detection of mini black holes at the LHC [26]. Thus, in gravity's rainbow, the last stages of the evaporation of a black hole is very different from general relativity. It has been explicitly demonstrated that such a remnant also occurs for black rings [27]. In fact, it has been proposed that such a black remnant will form for all black objects in the gravity's rainbow [28]. This has also been explicitly demonstrated for Kerr, Kerr-Newman-dS, charged-AdS, higher dimensional Kerr-AdS black holes and a black Saturn [28]. The black holes in the Gauss-Bonnet gravity have also been studied using gravity's rainbow [29]. It was demonstrate that even though the thermodynamics of the black holes get modified in the Gauss-Bonnet gravity's rainbow, the first law of thermodynamics still holds for this modified thermodynamics. As the Gauss-Bonnet gravity generalizes to Lovelock gravity, it is interesting to analyze the black hole solutions in Lovelock gravity using gravity's rainbow. Thus, in this paper, we will analyze certain aspects of black holes in Lovelock gravity using gravity's rainbow. Furthermore, string theory is one of the motivations for studding gravity's rainbow, and the low energy effective action for the string theory produces Lovelock gravity [30][31][32][33]. It means, in the context of string theory, higher curvature terms of Lovelock Lagrangian are string corrections on gravity side at the classical level. Thus, theoretically, it seems natural to study the effect of higher curvature terms in the context of gravity's rainbow. We also extend the investigation of Lovelock gravity's rainbow to the case of coupling with linear (and nonlinear) electromagnetic field. The Lovelock gravity coupled to different classes of the electromagnetic fields has also been studied [34][35][36][37][38][39][40][41][42][43][44][45][46]. It has been observed that this can have interesting phenomenological consequences [47]. Motivated by these subjects, we will analyze the black holes in Lovelock gravity's rainbow coupled to the Maxwell field, and then extend our investigations to the case of nonlinear electrodynamics.
We have organized this paper as follows: First, in Sec. II, we give a brief review of the proper field equations of Lovelock gravity in the presence of linear electromagnetic field, and then present a new class of charged black hole solutions in gravity's rainbow formalism. After that, in subsection II A, we study the effect of rainbow functions on thermodynamic quantities and show that the first law of thermodynamics still holds in the context of this formalism. Also, in subsection II B, we perform a thermal stability analysis for these black holes in the canonical ensemble. Next, in Sec. III, we extend our consideration to the case of Born-Infeld nonlinear electrodynamics (BI) and present a new class of charged black holes in Lovelock-BI gravity's rainbow, and then investigate their thermodynamic properties. In Sec. IV, in order to complete our discussion of thermodynamic properties, we study the critical behavior of black holes in the extended phase space by regarding the cosmological constant as a thermodynamical pressure. We finish our paper with some concluding remarks.

II. LOVELOCK-MAXWELL GRAVITY'S RAINBOW
Lovelock theory of gravity is one of the standard extensions of general relativity in higher dimensional spacetimes. One can consider the two first terms (zero and first term) of Lovelock Lagrangian to obtain a constant term related to the cosmological constant and also Ricci scalar, respectively. Higher curvature terms in the gravitational field equations appear as a result of adding the second and third terms (and also higher order terms) of Lovelock Lagrangian in the gravitational action. The theory obtained by considering first three terms of the Lovelock Lagrangian (and ignoring higher order curvature terms) is called Gauss-Bonnet gravity. Here, we want to expand our study to the third term (first four terms of Lovelock Lagrangian) in the presence of Maxwell electrodynamics; the so-called third order Lovelock-Maxwell gravity. The Lagrangian of the third order Lovelock-Maxwell gravity can be written as where L 0 = −2Λ and L 1 =R are, respectively, the cosmological constant and the Ricci scalar and α 0 = α 1 = 1. The second and third order Lovelock coefficients α 2 and α 3 indicate the strength of the second and third order curvature terms. Also, L 2 and L 3 are, respectively, the Gauss-Bonnet Lagrangian and the third order Lovelock term given as In addition, the last term of Eq. (1) is the Maxwell invariant F = F µν F µν , where F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic tensor related to the gauge potential A µ . Using the variational principle, the electromagnetic and gravitational field equations of third order Lovelock-Maxwell gravity can be obtained as where the cosmological constant term and Einstein tensor are, respectively, denoted by G µν and G µν are, respectively, the second and third orders Lovelock tensor given as Now, we consider the following line element to obtain the d-dimensional static topological black hole solutions in the context of gravity's rainbow where Using Eq. (5) the gauge potential is obtained as follows where q is an integration constant which is related to the electric charge. It is easy to show that the nonzero components of the electromagnetic tensor are Here, to have a practical solutions for gravitational field equations, we consider a special class, in which α 2 and α 3 are related to each other such as . The gravitational filed equations with this condition have one real and two complex (conjugate) solutions for the metric function Ψ(r) . The real metric function Ψ(r) has the following form where the parameter m is an integration constant which is related to the finite mass. It is obvious that the obtained solutions reduce to the Lovelock-Maxwell solutions with condition f (ε) = g(ε) = 1 (or correspondingly ε −→ 0). In order to investigate the possible curvature singularity, we can calculate the Kretschmann scalar. The Kretschmann scalar has the following form Inserting Eq. (12) into Eq. (13), one finds an essential singularity at the origin and for r = 0 all curvature invariants are nonsingular. Equation (13) confirms that the strengths of singularity and other finite values of curvature can be drastically affected by rainbow functions. Numerical calculations show that, depending on the values of parameters, the metric function has two real positive roots, one extreme root or it may be positive definite. Hence, obtained solutions can be covered with an event horizon and the solutions may be interpreted as the black holes with two horizons (Cauchy and event), extreme black holes or naked singularity. The asymptotical behavior of the solution (i.e. r → ∞) is obtained as follows and therefore these solutions are asymptotically AdS In order to study the effects of energy dependency, we have to use a specific functional form of the rainbow functions. There are different phenomenologically motivations for considering different forms of the rainbow functions. A specific form of the rainbow functions has been motivated from results obtained in the loop quantum gravity and noncommutative geometry [21,22]. The hard spectra from gamma-ray burster has also been used to motivate the construction of a different form of rainbow functions [48]. Another interesting form of rainbow functions is based on the modified dispersion relations, such that the constancy of velocity of the light is not violated [49], and it can be explicitly written as where λ is an arbitrary parameter. Hereafter, we focus on Eq. (15) as a typical form of rainbow functions. It may be noted that depending on the energy ratio ε, one can obtain different values for the rainbow functions (see table I). So, in this paper, we relax dynamic determination of these functions and analyze the behavior of the system from this set of phenomenologically motivate rainbow functions (motivated by the modified dispersion relation with constant velocity of light).

A. Thermodynamics of black holes in Lovelock-Maxwell gravity's rainbow
This section is devoted to calculation of conserved and thermodynamic quantities, and examination of the first law of thermodynamics. At first, we use the definition of surface gravity to calculate the temperature   where after simplification, we obtain Regarding Eq. (17), we find that black hole temperature is modified due to the presence of rainbow functions g(ε) and f (ε). Entropy is an extensive quantity corresponds to the temperature as an intensive quantity. In higher derivative gravity the area law of entropy is not valid, generally. Therefore, we calculate it by Wald method which is valid in higher derivative gravity such as Lovelock gravity. It is shown in third order Lovelock gravity the entropy can be written as whereR µνγδ ,R µν andR are, respectively, the Riemann tensor, the Ricci tensor and the Ricci scalar of the induced metricg µν on (d−2)-dimensional horizon. Calculation shows that the modified entropy of third order Lovelock gravity is obtained as follows in which shows that the area law is valid only for the Ricci flat solutions (k = 0). In order to establish the first law of thermodynamics for charged black holes, we have to calculate quantities related to the electrodynamics. First, in order to obtain its extensive quantity (i.e. the electric charge), we calculate the flux of the electromagnetic field at infinity. The electric charge is obtained as Now, we should calculate the electric potential of black hole solutions. Electric potential is an intensive quantity corresponds the electric charge, in which measured at the horizon with respect to infinity as a reference The finite mass of the black hole can be obtained by using different methods that all of them have the same result. Using the behavior of the metric at large r (for asymptotically flat black hole) or counterterm method (for asymptotically AdS black hole), it is easy to show that the finite mass of the black hole is Now, we rewrite the finite mass M as a function of the extended quantities (the entropy and electric charge). Straightforward calculations show that where r + = r + (S) and Now, we use the first law of thermodynamics to define temperature and electric potential as the intensive parameters conjugate to the entropy and electric charge Using a complete numerical analysis, it is easy to show that Eqs. (24) and (25) coincide with Eqs. (17) and (21), respectively, and therefore we deduce that all intensive and corresponding extensive parameters satisfy the first law of thermodynamics with the following form B. Thermal stability of black holes in Lovelock-Maxwell gravity's rainbow Here, we investigate local stability of the black hole solutions in Lovelock-Maxwell gravity's rainbow by using determinant of Hessian matrix. The local stability requires that the behavior of energy (i.e. M (S, Q)) be a convex function of its extensive quantities, S and Q. In canonical ensemble the electric charge is a fixed variable and so determinant of Hessian matrix is a function of the heat capacity as follows For allowed physical black hole parameters (such as positive temperature T and finite mass M ), the positivity of the heat capacity ensures the local stability. Analytical calculations of the heat capacity (or determinant of the Hessian matrix) are too large and, therefore, we leave out the analytical results for reasons of economy and brevity, and would rather offer a more practical discussion with some figures. Numerical calculations for thermal analysis show that for a black hole with definite mass M , there always is a lower critical value for the radius of event horizon, r +c . The black hole temperature is always positive for r + > r +c . Now, we can discuss the heat capacity of black holes. Taking into account the heat capacity, we find that there is a critical value α c for the Lovelock parameter. In the cases of allowable temperature, i.e. r + > r +c , black holes are thermally stable for large enough values of α, i.e. α > α c (see Fig. 1). In addition, we can find that for α < α c , there is an extreme radius for the event horizon (r +ext ) where in regions r +c < r + < r +ext the heat capacity is positive and therefore black holes are stable. Also, there is an upper limit for the radius of event horizon (r +u ) where in regions r +ext < r + < r +u black holes are unstable. Finally, for regions r + > r +u , we find a stable behavior for black holes (see Figs. 1-3 for more details).

III. LOVELOCK-BI GRAVITY'S RAINBOW
In this section, we generalize our consideration to the case of nonlinear electrodynamics (NED). We choose the Born-Infeld field (BI field) of NED to obtain the Lovelock-BI gravity's rainbow black holes. To construct this theory one can replace Maxwell Lagrangian (−F ) with the BI Lagrangian L(F ), which is a function of Maxwell invariant, F . The BI Lagrangian is where β is BI parameter with the dimension of (mass) 2 in natural units and as one expects, in the limit β → ∞, L(F ) reduces to the Maxwell one (−F ). Varying the following total Lagrangian with respect to the metric tensor g µν and gauge potential A µ , one can obtain the gravitational and electromagnetic field equations.
where α i 's and L i 's are defined before. Calculations show that the proper field equations are where L F = dL(F ) dF , and G (i) µν 's are defined in section II. In order to study Lovelock-BI gravity's rainbow, we use the defined rainbow metric (Eq. (8)), and find the gauge potential and the metric function as where in which H(η) is a hypergeometric function, and in the limit β → ∞ the obtained gauge potential and metric function reduce to those of the Maxwell gauge potential and Lovelock-Maxwell metric function (Eqs. (10) and (12)), respectively. The Kretschmann scalar diverges only at r = 0 for the new solutions, and therefore, there is an essential singularity at the origin. It is notable that the mentioned singularity can be covered with an event horizon and thus one can interpret the singularity as a black hole. It was shown that for nonlinearly charged black holes, depending on the values of α and β, the metric function may have two different behaviors; Reissner-Nordström like and Schwarzschildlike (for more details see [50,51]). In addition, for large values of r, these solutions reduce to previous Maxwell case, and therefore, the nonlinearity parameter does not change the asymptotical behavior. In order to investigate the effect of rainbow functions on the metric, we plot Fig. 4. According to this figure, one finds that both f (ε) and g(ε) not only affect the location of event horizon, but also their values characterize the type of horizon. In other words, decreasing rainbow functions leads to increasing the event horizon. In addition, regarding different values of rainbow functions, one can obtain timelike singularity with two horizons, an extreme horizon or naked singularity. Also there is a minimum value of rainbow function, in which for its lower values, singularity will be spacelike with a non-extreme event horizon.

A. Thermodynamics of black holes in Lovelock-BI gravity's rainbow
Here, we are going to check the first law of thermodynamics for the black hole solutions in Lovelock-BI gravity's rainbow. At the first step, we calculate the conserved and thermodynamic quantities. Surface gravity interpretation leads to the following Hawking temperature According to Eq. (34), we find that rainbow functions affect the Hawking temperature of the black holes. In addition, based on Fig. 5, one finds that changing rainbow functions affect the value of temperature and also increasing g(ε) (or decreasing f (ε)) leads to increasing the root of T . It is notable that the root of temperature is a bond point for the radius of event horizon, in which there is no physical black hole with horizon radius smaller than such bond point.
Entropy of the black hole can be calculated as the previous relationship which leads to the same relation as Eq. (19). The electric potential and charge of the black hole solutions are obtained as Since the finite mass of the solutions is the same as that of Sec. II, we can find the black hole mass as a function of the entropy and electric charge as where . Now, we are in a position to check the first law of thermodynamics. By computing ∂M/∂r + , ∂S/∂r + and ∂Q/∂r + and using chain rule, one can show that the following quantities are the same as the temperature and electric potential given in Eqs. (34) and (35), respectively. So, the obtained thermodynamic and conserved quantities satisfy the first law of thermodynamics, dM = T dS + ΦdQ.

B. Thermal stability of black holes in Lovelock-BI gravity's rainbow
Taking into account thermodynamic quantities, we are in a position to investigate the local stability of black hole solutions in Lovelock-BI gravity's rainbow. As we have mentioned before in Sec. II B, the local stability requires that the behavior of energy (i.e. M (S, Q)) be a convex function of its extensive quantities S and Q, and the positivity of the heat capacity guarantees the local stability in the canonical ensemble. Here, we continue our discussion of thermal stability with numerical analysis again. The only difference between this case and the case of Lovelock-Maxwell gravity's rainbow is adding the nonlinearity. In this case, the nonlinearity parameter affects the value of r +c . Calculations show that regardless of the values of β, there is always an unstable phase of black hole solutions for α < α c . So as before, there are two limits for the small enough values of α, in which black holes are unstable for r +ext < r + < r +u . Also, Numerical calculations for thermal analysis show that the black hole temperature is always positive for r + > r +c depending on the parameter β. Therefore, stability condition for α < α c shows that the heat capacity is positive for r +c < r + < r +ext and thus black holes are thermally stable in this region. Furthermore, for r + > r +u , we find a stable phase for black holes again. Finally, it is worthwhile to mention that for r + > r +c and α > α c , calculations show a thermally stable phase for the black hole solutions of Lovelock-BI gravity's rainbow (see Fig. 6 for more details).
For the sake of completeness, we plot Figs. 7 and 8 to investigate the effects of rainbow functions on thermal stability of the solutions. Left panels of these figures indicate the root of heat capacity which is coincident with the place of vanishing temperature. In addition, these figures show that there is a minimum value for the rainbow functions (f min (ε) and g min (ε)) which separates two different behaviors. In other words, black holes are thermally stable for f (ε) < f min (ε) (g(ε) < g min (ε)). Otherwise for f (ε) > f min (ε) (g(ε) > g min (ε)) there are two divergences for the heat capacity which indicate a phase transition. It is notable that the values of f min (ε) and g min (ε) can be numerically calculated as functions of other parameters.

IV. PV CRITICALLY OF LOVELOCK GRAVITY'S RAINBOW
In order to study the critical behavior of black holes in the extended phase space, we treat the cosmological constant as a thermodynamics pressure. In other words, we do not work in a fixed AdS or dS background. In fact, the cosmological constant is not a constant anymore, but a variable. According to the AdS/CFT correspondence, of interesting case is asymptotically AdS black holes, and in this section, we consider AdS black holes with spherically topological horizons (there is not any phase transition for k = 0 and k = −1 in this gravity model). If we treat the negative cosmological constant proportional to thermodynamic pressure, its conjugate quantity will be thermodynamic volume. This new assumption has considered by many authors in recent years (see [52][53][54][55][56][57][58][59][60][61][62][63][64][65] and references therein). Already, using scaling argument, it has shown that the Smarr relation is consistent with first law of thermodynamics by assuming the cosmological constant as a thermodynamic variable [46,61,[66][67][68]. Generalization of extended phase space thermodynamics to higher derivative gravity theories and investigation of triple points, reentrant phase transitions and equal area law have been studied before [43,44,[69][70][71][72][73][74][75][76][77][78][79]. Here, we are going to study the phase transition and critical behavior of black hole solutions in a more complicated gravitational background, i.e. for our interest, Lovelock gravity's rainbow. In the geometric units G N = = c = k B = 1, one can identify the cosmological constant with the pressure as It was seen that by considering the cosmological constant as thermodynamic pressure, the black hole mass M can be explained as enthalpy rather than internal energy of the system, i.e. M = H [66][67][68], and so the Gibbs free energy will be in the form of G = M − T S. Using the chain rule, thermodynamic volume of black hole is obtained as where X i denotes all other extensive quantities and w d−2 = V d−2 | k=1 is the volume of (d − 2)−dimensional unit sphere (see Eq. 9 with k = 1). The critical point in an isothermal P − V diagram is an inflection point, and therefore, it can be obtained from the following equations Investigation of the critical behavior is not possible, analytically, and so we have to use numerical analysis. In order to study the mentioned critical behavior, we present various tables and figures. It is worthwhile to mention that the characteristic swallow-tail form of G − T diagrams, inflection point of isothermal P − V plots and also subcritical isobar of T − V diagrams guarantee the existence of the phase transition. The associated P − V , T − V and G − T diagrams for both types of the obtained black hole solutions (Lovelock-Maxwell and Lovelock-BI) are displayed and the related critical points are obtained in various tables. The effects of higher derivative Lovelock gravity, nonlinearity parameter of BI theory and rainbow functions are evident in presented tables and figures.  Table II: Lovelock-Maxwell solutions: k = 1, q = 1, f (ε) = 0.9, g(ε) = 0.9 and d = 7.   Table III: Lovelock-Maxwell solutions: k = 1, q = 1, f (ε) = 0.9, g(ε) = 0.9 and α = 0.5.  Table IV: Lovelock-Maxwell solutions: k = 1, q = 1, g(ε) = 0.9, d = 7 and α = 0.5.  Table VI: Lovelock-BI solutions: k = 1, q = 1, f (ε) = 0.9, g(ε) = 0.9, d = 7 and α = 0.5.
In Figs. 9-14 we show the critical behavior of the system in both Lovelock-Maxwell and Lovelock-BI gravity's rainbow. Also, we present five tables to investigate the effects of various parameters on the critical point. According to these figures and related tables we find that the Lovelock parameter (α), nonlinearity parameter (β), dimensionality (d) and rainbow functions affect the critical point, significantly. In addition, variation of these parameters can change the near universal ratio PcVc Tc , considerably. Regarding Figs. 13 and 14 with tables IV and V, we can find the effects of changing rainbow functions with more details. Based on Fig. 13 and table V, we find that the critical values of pressure and temperature (volume) increase (decreases) by increasing g(ε). While Fig. 14 and table IV show that increasing f (ε) leads to decreasing all critical values. Looking at the relation PcVc Tc , we find that it is an increasing function of f (ε), as opposed to increasing g(ε) where such ratio is a decreasing function. Therefore, one can adjust the values of these parameters to cancel the effects of all free parameters and obtain the universal ratio.

V. CONCLUSION
In this paper, we have analyzed the thermodynamics of black holes in Lovelock gravity's rainbow. We have obtained black hole solutions for both Lovelock-Maxwell and Lovelock-BI gravity's rainbow theories with different horizon topologies. We have also found that the asymptotical behavior of the solutions may be (a)dS or flat with an effective energy dependent cosmological constant. We have calculated conserved and thermodynamic quantities in the energy dependent background and found that although rainbow functions may change these quantities, the first law of thermodynamics is valid. In addition, we have examined thermal stability of the solutions in the context of canonical ensemble. We have shown depending on the chosen parameters, there is a critical horizon radius (r +u ), in which for r + > r +u the black holes are thermally stable.
At last, we have investigated the critical behavior and phase transition in the extended phase space by considering the proportionality between Λ and pressure. We have studied the effects of α, β, d and rainbow functions on the critical values and found that they can change critical point, significantly. In other words, we have shown that in addition to Lovelock and nonlinearity parameters, rainbow functions affect thermodynamic nature of the black hole. We have found that depending on the values of rainbow functions, a phase transition can occur. While one could adjust rainbow functions to obtain thermally stable black holes. In brief, we should note that rainbow functions affect various aspects of a black hole such as: strength of curvature and singularity, place and type of event horizon, asymptotical behavior, thermodynamic quantities, thermal stability and existence of phase transition.
It may be noted that various interesting systems have been analyzed using Lovelock gravity. The quasinormal modes of black holes in Lovelock gravity has been studied [80]. In this study, the WKB method was applied for analyzing the quasinormal frequencies Lovelock gravity. It would be interesting to analyze the quasinormal modes of black holes in Lovelock gravity's rainbow. It is expected that the choice of the rainbow functions will effect the behavior of these quasinormal modes. A scalar-tensor version of Lovelock theory with a non trivial higher order galileon term involving the coupling of the Lovelock two tensor with derivatives of the scalar galileon field has been construed [81]. It would be interesting to construct this theory using Lovelock gravity's rainbow, and then use it for analyzing black hole solutions. It has also been demonstrated that a black remnant forms in gravity's rainbow [25], and this has important phenomenological consequences for the detection of mini black holes at the LHC [26]. It would be interesting to investigate the formation of such black remnants in Lovelock gravity's rainbow. In fact, it has been demonstrated that such remnants will form for all black objects in usual gravity's rainbow [28]. It will be interesting to investigate if such a result holds for all black holes in Lovelock gravity's rainbow.