On generalized quasi-topological cubic-quartic gravity: thermodynamics and holography

We investigate the thermodynamic behaviour of asymptotically anti de Sitter black holes in generalized quasi-topological gravity containing terms both cubic and quartic in the curvature. We investigate the general conditions required for physical phase transitions and critical behaviour in any dimension and then consider in detail specific properties in spacetime dimensions 4, 5, and 6. We find for spherical black holes that there are respectively at most two and three physical critical points in five and six dimensions. For hyperbolic black holes we find the occurrence of Van der Waals phase transitions in four dimensions and reverse Van der Waals phase transitions in dimensions greater than 4 if both cubic and quartic curvature terms are present. We also observe the occurrence of phase transitions in for fixed chemical potential. We consider some applications of our work in the dual CFT, investigating how the ratio of viscosity to entropy is modified by inclusion of these higher curvature terms. We conclude that the presence of the quartic curvature term results in a violation of the KSS bound in five dimensions, but not in other dimensions.


Introduction
In recent years there has been an increased appreciation of the role that higher-curvature theories of gravity play in our understanding of several areas of physics, including supergravity and string theory, holography, cosmology, and black holes. These studies are motivated both by a desire to understand the ultraviolet behaviour of gravity and by a realization that terms non-linear in the curvature generically appear in perturbative calculations, particularly in string theory [1]. Furthermore, an analysis of this class of theories provides us with new insights into general relativity (or Einstein gravity) and may even provide new empirical tests of gravitational physics [2][3][4][5][6][7][8][9][10][11].

JHEP07(2019)012
Originally proposed nearly a century ago [12,13], a revival of interest in these theories in the theoretical physics community occurred once significant effort began to be expended on constructing a quantized theory of gravity. Adding terms quadratic in the curvature to the Einstein-Hilbert action were found to yield a power-counting renormalizable theory [14], and later a Gauss-Bonnet term was found to appear in the low energy effective action of string theory [15]. More recently it has been shown from a variety of perspectives [16][17][18] that a proper investigation of dual theories beyond large N in the context of the AdS/CFT correspondence conjecture [19,20] entails inclusion of higher curvature terms.
A key challenge presented by a generic higher-curvature theory of gravity is that the equations of motion are greater than second order in the derivatives, leading to a number of pathological properties such as the appearance of ghost degrees of freedom and other instabilities. However there exist a few classes of theories in which such pathologies are ameliorated and in certain cases are absent. The best known example is the Lovelock class of gravitational theories [21]. This class yields second order equations of motion in arbitrary dimensions, with the Einstein-Hilbert term being one of several terms that constitute Lovelock theory in a given dimension. In this sense Einstein gravity can be regarded as a special case of Lovelock gravity. Since this class of theories is ghost-free [15] they are viable candidates for generalizations of Einstein gravity in higher dimensions d ≥ 2k + 1 for a Lovelock theory that is kth order in the curvature. For dimensions d < 2k + 1 such terms play no role in the equations of motion. Hence only k = 1 Einstein gravity has non-trivial field equations and so one must look beyond Lovelock gravity to obtain interesting higher curvature theories with implications in (3 + 1) dimensions.
In the past several years considerable progress has been made along these lines. A broader class of quasi-topological gravity theories [16,22] have been constructed that retain many of the nice properties of Lovelock gravity under certain symmetry restrictions. They are non-trivial in any dimension d ≥ 5 regardless of the order in the curvature. Cubiccurvature quasi-topological gravity, for example exists in d ≥ 5 whereas cubic Lovelock gravity exists in d ≥ 7. Furthermore, their field equations, while generally greater than second-order, become second order under the imposition of spherical symmetry. Moreover, the linearized of the equations of motion of quasi-topological gravity coincide with those of Einstein gravity on maximally symmetric spacetime backgrounds up to an overall prefactor [23], ensuring that negative energy excitations do not propagate to asymptotic regions of constant curvature.
Even more recently a more general class of higher-curvature gravity theories have been discovered that are of considerable interest both holographically and phenomenologically. This is because they are free of ghosts on constant curvature backgrounds, solutions of their field equations yield a metric that depends on a single metric function in the spherically symmetric case, and they are dynamically non-trivial even in four dimensions. First obtained in (3 + 1) dimensions for cubic curvature [24] (a theory known as Einsteinian cubic gravity or ECG), they were found to have generalizations to any dimension [25] and to quartic powers in the curvature [26]. Generalizations to any power in the curvature exist [27] but have not been explicitly constructed.

JHEP07(2019)012
This class of theories -referred to as Generalized Quasitopological Gravity (GQG) -has several remarkable features. The constraint that spherically symmetric solutions depend on only a single metric function is found to also eliminate ghosts upon linearization of the theory on a constant curvature background [25]. Very recently it has been shown that they have a well-posed initial value problem for cosmological solutions and the potential to provide a late-time cosmology arbitrarily close to the Λ-Cold dark matter scenario whilst having a purely geometrical inflationary period in the early universe with a graceful exit [28][29][30]. While the field equations can be solved exactly in certain special cases [31], it is possible to analytically investigate the thermodynamics of black holes even in the generic case where analytic solutions are not available [32,33]. Charged black branes have an interesting phase structure that is absent for both their Lovelock and quasi-topological black brane counterparts [34]. The Kovtun-Son-Starinets bound on the ratio of entropy density to shear viscosity always holds [35], and small asymptotically flat black hole solutions were found to be stable [27], which may have implications for the information loss problem. Further holographic applications of this class of theories have been carried out, with discussions of the a-theorem, the universal stress-tensor two-point function and a universal relation for central charges [36][37][38]; more recently an extensive investigation of the holographic properties of the cubic case without massive modes was carried out, including holographic central charges, energy flux, Renyi entropies, and shear viscosity to entropy ratio [39]. Other recent work has shown that the shadows of GQG black holes have potentially interesting phenomenological signatures [9,10].
Thermodynamics of GQG black holes has yet to be fully explored. Previous studies for asymptotically flat solutions and AdS black holes have appeared in restricted contexts [25,27,34,40], but a full study combining both cubic [25] and quartic GQG [26] has yet to be carried out. The purpose of this paper is to conduct such a study for both spherical and hyperbolic charged black holes.
Our investigation will be carried out in the context of black hole chemistry, in which the cosmological constant is taken to be a thermodynamic variable [41,42] that is interpreted as pressure in the first law of black hole mechanics [43,44]. An extensive amount of work over the past six years has been carried out in this subject [45] and has indicated that black holes can exhibit a broad range of phase behaviour that has been observed in other areas of physics. Examples include triple points [46], re-entrant phase transitions [47], polymer-like behaviour [48], and even superfluid-like phase transitions [49][50][51], as well as a deep analogy between charged anti-de Sitter black holes and Van der Waals fluids [52]. Higher-curvature gravity theories have, using this approach, likewise been seen to have a very rich thermodynamic structure [33,[48][49][50], with even more results surveyed in a recent review [45].
Our paper is organized as follows. In section 2 we present charged static, spherically symmetric AdS black holes in the cubic-quartic GQG theory. This includes an asymptotic solution, a near horizon solution and then their match in the form of a numerical solution. In section 3 we investigate the thermodynamic properties of charged black holes in cubic-quartic generalized quasi-topological gravity by applying the black hole chemistry formalism. In section 4 we classify the phase structure and critical points for these black JHEP07(2019)012 holes by considering the perspective of black hole chemistry, working in the fixed charge ensemble. In section 5 we consider the thermodynamics for fixed potential ensemble. We will analyze the four dimensional case in detail, and then present relevant results in higher dimensions. In section 6 we present some results of holographic hydrodynamics to understand these theories in the context of AdS/CFT correspondence. We summarize our work in section 7 and present some directions for further research.

Charged black hole solutions in cubic-quartic GQG
We begin by setting up charged static, spherically symmetric AdS black holes obtained from the equations of motion that follow from a combination of cubic and quartic terms of generalized quasi-topological gravity (GQG).

Construction of equations of motion
Consider the class of static radially symmetric metrics with radial coordinate r and time coordinate t. Choosing coordinates so that the radius of a (d − 2) sphere behaves as r d−2 , higher-curvature gravity theories up to quartic order are obtained by applying the condition g tt g rr = −1 in order to get a single metric function. They result in both Lovelock and quasi-topological curvature terms as well as GQG curvature terms.
Concentrating only on the properties of GQG theory, we put aside the Lovelock and quasi-topological terms (for a discussion see e.g. [62,66]) and consider Einstein gravity accompanied by cubic and quartic generalized quasi-topological terms, with minimal coupling to an Abelian gauge field. The action 1 in d dimensional spacetime reads [26] where the cosmological constant Λ = − (d−1)(d−2) and the quartic generalized quasi-topological term [26] is given at appendix A. The rescaled cubic couplingμ, and quartic couplingλ, are given bŷ where µ and λ are arbitrary coupling constants, and rescaling is done to simplify the field equations.

JHEP07(2019)012
As per our requirements for radially symmetric metrics, we employ the following ansatz and we find that the field equations of GQG yield N (r) = constant [25]; we shall set N (r) = 1 for simplicity. 2 Here dΣ 2 (d−2),k describes the (d−2)-dimensional line element of the transverse space, where k = +1, 0, −1 stand for spherical, flat and hyperbolic geometries of a surface of constant scalar curvature. As an investigation of the k = 0 case has previously been carried out [34], we shall in the sequel consider only non-planar black holes.
For a maximally symmetric space, the metric (2.4) becomes, where is the AdS radius that is related to the cosmological constant. The quantity f ∞ = lim r→∞ f (r) 2 /r 2 and is obtained by solving following polynomial equation which is independent of the choice of k in the transverse section. While at least one coupling is non-zero, the higher curvature terms drive away f ∞ from unity. Since we require the same asymptotics as AdS space we only pick positive real solutions of the above polynomial. The effective radius of the AdS space is given by eff = / √ f ∞ . In fact it turns out the negative of the derivative of eq. (2.6) with respect to f ∞ yields the prefactor of the linearized equations of motion [25], and to prevent the appearance of ghosts in the particle spectrum we require P (f ∞ ) > 0. For charged black holes, we include a Maxwell field F ab = ∂ a A b − ∂ b A a , with electromagnetic one form defined as A = qE(r)dt, (2.8) where, by inserting the above expression into the Maxwell equation, we get for the electric field. The specific choice of prefactor makes for greater simplification later on in the field equations; we choose the constant term in the potential to be zero. 2 We note that the choice N = 1/ √ f∞ has been used [16] to normalize the speed of light on the boundary or to get c = 1 in the dual CFT. Here we shall set N = 1, and note that by time reparametrization of the metric we can obtain c = 1 on the boundary if desired.

JHEP07(2019)012
The field equation for the action (2.1) yields the relation where m is a constant of integration and where [25] 12) are respectively generated by the cubic and quartic generalized quasi-topological terms.
The parameter m has scaling dimension [length] d−3 and we will see that appears in the formula for the mass of black hole. While exact solutions to the field equation seem hard to find (except in special cases [31]), studying the far and near horizon behaviour of the metric perturbatively is still feasible, and permits us to analytically obtain the thermodynamic quantities associated with black hole solutions. Specifically, we shall utilize information from the near horizon expansion to describe the thermodynamics of the black holes.

Degenerate vacuum solutions
Before discussing solutions with nonzero mass and charge, we first consider solutions with degenerate vacua. These are analogous to the Lovelock-Unique-Vacuum (LUV) solutions [78][79][80][81], and occur when the field equation (2.10) has solutions of the form (2.5) with multiple degenerate solutions for L = / √ f ∞ for q = m = 0, or alternatively when (2.6) has multiple degenerate solutions.
It is straightforward to show that there is a one-parameter family of doubly-degenerate solutions to (2.10) for q = m = 0, given by

JHEP07(2019)012
valid for any f ∞ and . There is also a triply-degenerate solution 14) with f ∞ = 2. Note that for d = 8 and d = 6 these solutions do not exist. There are no fully degenerate solutions of the LUV type to (2.10). This situation could presumably be altered if we were to add the quadratic Lovelock term (the Gauss-Bonnet term) to (2.1), but we shall not pursue this here.

Far region solution
In the asymptotic limit, the form of the metric function is Inserting the above expansion into eq. (2.10) and requiring that it be satisfied at each order in a 1/r expansion yields where P (f ∞ ) was introduced in (2.7). We have presented the six leading terms, and displayed the schematic structure of the next order corrections to f (r) asymp . It is obvious that for µ → 0 and λ → 0 that P (f ∞ ) → 1 and so as expected from the solution in Einstein gravity.

JHEP07(2019)012
Degenerate (multiple) solutions of (2.6) will yield weaker falloff behaviour for nonzero m and q, as is the case with LUV-type solutions [81]. Setting m = q = 0, f = f ADS turns out to be an exact solution to the cubic theory [31] and is also an exact solution to the field equations of the cubic-quartic theory (2.1) provided (2.6) holds. The thermodynamic behaviour of this class of solutions is uninteresting (essentially it is the same as that of the BTZ black hole [82]). We shall briefly consider this situation in section 4.1, where we shall see it is connected with a maximal pressure.
To find the asymptotic behaviour in the degenerate case, and for simplicity, we consider the cubic and quartic parts of the action separately. In both cases, only for doubly degenerate vacuum solutions (referred to as the critical state) in which P (f ∞ ) in (2.7) vanishes do we get nonvanishing coupling. For higher orders of degeneracy (where all derivatives of eq. (2.6) are zero up to a specific order) we find vanishing coupling.
For the cubic case, where the associated critical coupling is given in (4.5) and f ∞ = 3/2, higher derivatives of P (f ∞ ) at this point are nonzero up to the third derivative of (2.6). Studying the asymptotic behaviour for degenerate vacua can be carried out by inserting into the equation of motion (2.10) with vanishing quartic coupling for now. For the parameters at the critical state we get a differential equation for (r). To find an asymtotic solution to this equation we write where (r) is the correction to the asymptotic AdS metric, and A and x are determined upon insertion into (2.10). Setting q = 0 (for simplicity), we find showing (as before) that in six dimensions there is no non-vanishing solution. For other dimensions there are two values for A allowed (both signs of the square root). Our result for x is compatible with the result given in [83] in the case of double degeneracy. Finding subleading terms requires including more terms in the initial expansion of (r). For the quartic case the critical state only exists at the quartic coupling given in (4.9) and f ∞ = 4/3, which yields two degenerate vacuum states while P (f ∞ ) vanishes and up to forth derivative of eq. (2.6) are nonzero at this critical value. Repeating a similar procedure for the cubic part, at leading order for large r and A = 0 we obtain where for d = 8 there is no nonzero solution. For d < 8 there are two solutions for A; for d > 8 there is no real value for A unless m < 0. We obtain the same falloff for the radial coordinate as in the cubic case, since it governs the value of the degeneracy.

JHEP07(2019)012
A more detailed study of 'excitations' of the degenerate cases (2.13) and (2.14) with nonzero m and q are somewhat more complicated than the Lovelock case since (2.10) is no longer a simple polynomial; we leave this for future investigation.
The homogeneous solution in the far region is found by inserting f (r) = f (r) asymp + f h (r) in eq. (2.10), where parameterizes the strength of these corrections. Substituting this expression in the field equation, we get an inhomogeneous second order differential equation for the function f h (r). At leading order in , assuming that µ = 0 and λ = 0 the homogenous part of the equation at large r becomes 3 and we note that it is independent of the value of k; for vanishing µ and λ it yields the well known AdS Reissner-Nordstrom (RN ) solution. For d = 6 there is an ambiguity in the above expression for µ = 0, λ = 0. For this particular case the applicable equation becomes 6 and the explicit form of γ 2 can be read from (2.24). The solution of (2.22) in the case of γ 2 > 0 is where I and K are the modified Bessel functions of the first and second kinds, and A and B are constants. In the limit of large r f h+ ∼ Ar 5/2 exp 2γr and so we must set A = 0 to ensure the AdS boundary conditions are satisfied. As a result no ghost excitations can propagate to infinity. We shall see shortly that the contribution of the second term can be dismissed. The solution for the particular case (2.24) can be obtained in a similar fashion using its corresponding value for γ 2 . Notice that k does not appear in the asymptotic solution and the numerator of γ 2 in conjunction with the positivity condition (2.7) ensures freedom from ghosts [25,26]. Indeed, the positivity of the numerator relation (2.23) gives the same no-ghost condition as in (2.7), and one only needs to check whether the denominator is positive as well.

JHEP07(2019)012
To have the correct asymptotics we must choose (µ, λ) and the mass parameter in such a way that γ 2 > 0. To see this, note that if γ 2 < 0 the homogenous solution asymptotically takes following form where J and Y are respectively Bessel functions of the first and second kind. In this situation, in any dimension the solution oscillates rapidly and its amplitude becomes larger than r 2 2 at large r. It therefore does not approach AdS asymptotically, and so we set C 1 = C 2 = 0 to get rid of this homogenous part of the solution. For the rest of our considerations, to avoid any oscillating behaviour near infinity we restrict the solutions to the constraint γ 2 > 0. Finally we note that the particular solution (2.16) polynomially decreases with 1/r, and is the dominant part of the total solution f (r) = f h+ + f (r) asymp for sufficiently large r; we therefore neglect the term f h+ in eq. (2.25) in the sequel.

Near horizon solution
To construct the solution near the horizon we consider the following expansion for the metric function, where T is the Hawking temperature of the black hole. It is found by imposing the regularity condition for the Euclidean sector of the complex manifold (under t → iτ ) and reads as Substituting the near horizon expansion of the metric function into the field equation and imposing that it holds at each order of (r−r + ) we obtain for the zeroth and first order terms which specify the formula for the mass parameter m and temperature T in terms of the horizon radius and coupling constants. We shall use these equations for our thermodynamic investigation later on. Continuing to higher order terms one is able to find all other series coefficients in terms of a 2 , which is a free parameter and its value is determined by using the boundary condition at infinity. We pause to comment that in a general non-linear theory of gravity a spherically symmetric metric depends on two functions, and the mass and the temperature are determined by two parameters. However it is a special feature of the generalized quasi-topological theories (as well as the Lovelock gravity theories) that both the mass and the temperature are determined in terms of one parameter: the horizon length, as (2.30) and (2.31) indicate. This is an important feature of these theories, since it allows us to analyze thermodynamic behaviour without full knowledge of the solution [26,32,33,35,40]. Now that we constructed the asymptotic and near horizon solutions, we next find a numerical solution that interpolates between them. For this purpose we define the rescaled metric function where g(r) → 1 as r → ∞. Here f ∞ is a positive real root of eq. (2.6). Choosing some specific values for the coupling constants, electric charge and horizon radius, while = 1, we find the associated values of mass parameter and temperature referring to (2.30) and (2.31).
To numerically solve the second order differential equation (2.10) we need to identify initial values for the metric function f and its first derivative. We use the value of f close to the horizon to set up the seed solution for g. We then fix the value of a 2 to desired order, using the shooting method such that the numerical solution for g approaches unity asymptotically.
As there are several branches of solutions, we select the one that tends to the Reissner-Nordstrom (RN ) solution in the µ → 0, λ → 0 limit; otherwise we get other solutions that are not physically interesting. Furthermore, because the differential equation is stiff, the JHEP07(2019)012 solution can be obtained only to a certain precision. For our choice of a 2 the asymptotic solution up to O(r −12 ) is precise to one part in 1,000 or better.
In order to exhibit the behaviour of the solution while varying the cubic and quartic couplings individually, we performed the computation for the cubic and quartic parts of the equation separately to see the impact of each of these terms individually. Figure 1 illustrates the numerical solution in four dimensions and shows that at fixed quartic coupling, increasing charge drives the horizon inward, whereas at fixed electric charge, larger values of |λ| displace the event horizon outward. 4 The right panel demonstrates the difference between the numerical result for the metric function and its corresponding asymptotic behaviour from (2.16); we see that these converge at enough large r. We also plotted the corresponding graphs for cubic gravity and find that similar behaviour takes place as charge and/or the cubic coupling are varied.
To see the behaviour of the metric function in four dimensions as r → 0, we expand the field equations in powers of r. In general we have the following expansion f (r) = r s (a 0 + ra 1 + r 2 a 2 + · · · ). (2.33) Considering only the quartic curvature part of (2.10), for small r the first term in the above expansion is the dominant contribution to the metric function. To find s, we use the same numerical procedure described before, and depict the behaviour of rf (r)/f (r) near r = 0. We find that in four dimensions s vanishes. However in higher dimensions depending on the choice of parameters s gains an non-integer negative value, except in six dimensions where it becomes positive for the choice of physical parameters as prescribed in the next section. A similar argument was given in [84] for the cubic part. Therefore in four dimensions the metric is regular at the origin. However, the Kretschmann scalar R abcd R abcd ∼ r −4 ; the spacetime is still singular, but its singularity is softer than its counterpart in Einstein gravity, in which R abcd R abcd ∼ r −6 .
We also find that the associated plots for the five-dimensional metric function (2.32) are similar, provided the parameters are chosen to satisfy a physical constraint that we will discuss in the next section.

Thermodynamic properties
Our aim is to study the effects of including cubic and quartic generalized quasi-topological terms on the known behaviour Einstein black hole thermodynamics in four and higher dimensions. We begin by investigating the first law and Smarr relation, where we apply the black hole chemistry formalism [45] by taking the cosmological constant, Λ and the couplings µ, λ as thermodynamic variables. We then study the physical constraints and find out whether at different domains in terms of the couplings and the charge they satisfy these constraints. We also elucidate the critical behaviour for these black holes. The red line shows the solution for Einstein gravity (for which the quartic coupling λ = 0) with zero electric charge. For the other curves we choose different values for these parameters. Top right: this panel depicts the difference between the numerical solution and its corresponding large-r analytic solution in quartic gravity -we see that convergence holds asymptotically. Bottom left: the plot presents the rescaled metric function (2.32) in four dimensions for cubic gravity (with λ = 0); the red line is again the solution for Einstein gravity with zero charge. Bottom right: the difference between the numerical solution and its corresponding large-r analytic solution in cubic gravity; we see again the asymptotic convergence. In all cases, we set = 1 and r + = 10. The mass parameter m is defined in (2.30).

First law and Smarr relation
As discussed in section 2.4, equations (2.30) and (2.31) provide the relations for obtaining the mass and temperature of the black holes without requiring knowledge of an exact solution. Since the explicit form for the temperature is complicated, we apply the second equation implicitly to verify the first law of thermodynamics is satisfied.
We utilize the Iyer-Wald formalism [85,86] to compute the entropy andε ab is the binormal to the horizon, which is normalized asε abε ab = −2. The induced metric on the horizon is γ ab and γ = detγ ab . From the action (2.1) we find for the entropy, where Σ (d−2),k is the volume of the submanifold with line element dΣ (d−2),k . For k = 1 this is the volume of the (d − 2)-dimensional sphere and finite, although when k = 0 and k = −1 one needs to perform some kind of identification to define finite volume.
Identifying the pressure as [45,87] the other thermodynamic quantities are and [88] the mass is It is straightforward to show these quantities satisfy the extended first law of black hole thermodynamics,

JHEP07(2019)012
where V is the thermodynamic volume conjugate to the pressure, and Ψ µ , Ψ λ are the respective thermodynamic conjugates to the couplings µ, λ. Furthermore, a scaling argument [44] applied to the various thermodynamic quantities above yields the Smarr formula which can be shown to hold for these quantities.
To investigate the critical behaviour of these black holes, an equation of state is required. This is obtained by substitution of 2 in (2.31) in terms of pressure. Hence where the different parameters are where v is the specific volume [89]. As in previous studies [33,34], we see from (3.9) that there is a non-linear dependence of the equation of state on the temperature. For future reference we choose the free parameters to be e, β 3 and α 4 ; these are independent of the choice of k.
The explicit form of the Gibbs free energy as where we pulled out an overall positive factor to simplify the expression; the explicit form of the other parameters is given in eq. (3.10). The equilibrium state is the one that minimizes the Gibbs free energy G for fixed temperature and pressure.

Physical constraints
We now explicate the constraints on the cubic and quartic couplings required for physical solutions. Generalized quasi-topological theories have the property that only the massless graviton propagates on constant curvature backgrounds provided the parameters are appropriately constrained. To ensure this, the effective Newton constant of gravity must have the same sign as that in Einstein gravity. This implies that the pre-factor in the linearized equations of motion about the AdS solution is positive [25], i.e., P (f ∞ ) > 0 with P (f ∞ ) defined in eq. (2.7) and the value of f ∞ is given by solution of eq. (2.6) which is positive in order to get an asymptotic AdS solution. The same relation occurs if we require γ 2 > 0 (see the discussion after (2.26)).
In terms of the rescaled parameters in eq. (3.10) and the pressure given in (3.4), the no-ghost constraint (2.7) becomes and we note that in the limit β 3 → 0, α 4 → 0 (or µ → 0, λ → 0) that we reach the Einstein branch of the theory. Disregarding solutions with γ 2 < 0 (since they are not asymptotically AdS) we have from (2.23) upon using eq. (3.10), where P (f ∞ ) is given in (3.12). It is well-known that in higher curvature gravity black hole entropy can be negative in some regions of parameter space, perhaps indicative of an occurrence of instability [90,91]. Imposing the requirement for positive black hole entropy yields (3.14) When temperature and specific volume are positive, in each dimension the values of couplings must be chosen to satisfy the above inequality. We search for the domains in parameter space where these conditions are valid for various charge and coupling constants in the next section.

Thermodynamics in the canonical ensemble
Equipped with the field equations and relevant thermodynamic relations, we consider first the fixed charge ensemble. We aim to investigate the phase structure and critical points for these black holes.

JHEP07(2019)012
The equation of state in terms of the rescaled parameters e, β 3 and α 4 introduced in eq. (3.10) is for any d ≥ 4; we note that only e 2 appears everywhere, so our results are valid for both positive and negative charge. Phase transitions occur if the equation of state demonstrates some oscillatory behaviour, with P (v) having at least one minimum and one maximum. This in turn depends on the signs of the coefficients of different powers of v, as these determine how many roots exist in the equations for critical volume and temperature. To get a critical point, the following equations must hold We also find that when µ = 0 in four and five dimensions λ must be negative, whereas in higher dimensions λ must be positive in order to get physical points that satisfy all physical constraints mentioned in the previous section. To get the critical volume and temperature in terms of charge and couplings, we solve equation in various dimensions. As the explicit form is lengthy, we do not present results explicitly; in practice it is easier to solve the equations parametrically for T and v c in terms of the other parameters in certain dimensions.

AdS 4 vacua and maximum pressure
We consider here the structure of the AdS 4 vacua of (2.1) in four dimensional spacetime, with curvature scale 1/ 2 eff = f ∞ / 2 . Setting the action length scale to = 1 (implying a fixed pressure of 3/(8π)), we analyze solutions to (2.6), considering the cubic and quartic couplings in distinction for simplicity.
Starting with only a non-zero cubic coupling µ, we have both µ < µ c and µ > 0; the first of these yields a negative real valued solution, and the second implies γ 2 < 0 in (2.23). Both regions are unphysical and we exclude them from further analysis. Note that from the linearized equations of Einstein gravity, the following relation holds between the effective Newton constant G eff and G. We see that for f 2 ∞ = − 4 4032µ , G eff → ∞, and inserting this value inside eq. (4.3) yields the critical coupling of µ c (noted previously [31]). At this point the discriminant of the cubic changes sign and there are distinct branches of solutions for values of µ < µ c and µ > µ c .
Repeating the same approach for higher dimensions, the overall behaviour of f ∞ given in figure 2 is similar to that in four dimensions. The critical limit is given by More generally we can reconsider the above discussion for arbitrary values of the pressure P . In four and five dimensions there is a maximum value for the pressure that results from the condition that the discriminant ∆ > 0, which is a constraint on parameter space JHEP07(2019)012 for physically acceptable solutions. In general we have where in four and five dimensions β 3 > 0 and is given in (3.10). For d = 6 the pressure is unbounded and for d ≥ 7 a similar procedure does not yield an upper bound for the pressure. Turning now to the quartic case, equation (2.6) becomes in d = 4 with = 1. The right graph in figure 2 indicates three possible real solutions to the above equation. The discriminant of (4.7) vanishes for λ = λ c = −81 6 /1024. Again, we note that G eff becomes infinity or equivalently P (f ∞ ) vanishes for f 3 ∞ = −3 6 /(16λ) and that (4.7) in turn implies λ = λ c . Similar to the previous case, for λ c < λ < 0, we have ∆ < 0 and there are two positive real solutions for f ∞ . Only the smaller of these has positive P (f ∞ ) and γ 2 . For λ < λ c , ∆ > 0 and there are no real solutions for f ∞ ; for λ > 0, although there is one real positive solution to (4.7), it implies γ 2 < 0 and is therefore physically inadmissible.
Requiring λ c < λ < 0 we conclude that there is a maximum value for the pressure given by where α 4 is given in (3.10) and Only for d = 4, 5 is α 4 > 0 and P max > 0, yielding an upper bound on the pressure; in higher dimensions there is no bound on the pressure.
In what follows, we concentrate on several specific dimensions and investigate the thermodynamic behaviour in some details. In the figures 3, 6, 11 given in the following sections we will adhere to the colour coding explained in table 1.

Critical behaviour in four dimensions
Our next task is to determine how many of these possible critical points are actually physical, and study their critical behaviour. We proceed by examining each value of d in succession.
Here we consider the effects of both cubic and quartic GQG in d = 4. The equation of state (3.9) becomes (4.10)

JHEP07(2019)012
It is obvious that for small v (i.e. for small black holes) that the term cubic in T (coming from quartic GQG) dominates. By taking different linear combinations of (4.2) it is possible to obtain an equation linear in T ; the resultant critical temperature and volume are then easily seen to satisfy the equations and which can be solved numerically for any choice of parameters for T c , v c . For simplicity, consider the behaviour of the critical temperature and volume, with only the quartic coupling active. Equations (4.11) and (4.12) become , and the critical volume v c satisfies We see that critical temperature is singular if the critical volume is such that the polynomial quartic in v 2 c in the denominator vanishes. An exception to this is if α 4 = 6400/729π 6 e 6 : the numerator in (4.13) also vanishes and t c remains finite. Note that is only occurs if k = 0, in accord with earlier work on black branes in GQG [34]. However the corresponding T c and v c become imaginary in the case of k = −1, so for hyperbolic black holes this singularity is absent.
For k = 1, if β 3 = 0 and α 4 = 6400/729π 6 e 6 , the resulting values of T c , v c correspond to a critical point with standard critical exponents (see (4.15) and the discussion following), and the phase transition in the vicinity of this point is a standard first order VdW transition similar to what is depicted in figure 4.
In studying the behaviour at this point, we note that the equations of state need to be solved for these specific parameter values, instead of using (4.13) and (4.14), since the latter becomes invalid if the denominator of t c vanishes. If both couplings are non-zero, again the denominator of T c in (4.11) is quartic in v 2 , and a similar procedure can be employed to write a formula for α 4 in terms of β 3 and e. Doing so, we find that for any values of the parameters the only solution is α 4 = 0, however in cubic gravity critical temperature does not have a singularity in four dimensions [84]. In other words, this particular occurrence of this apparent thermodynamic singularity is obtained only in Einstein-quartic GQG (for which the cubic coupling vanishes).  Due to the complexity of the equation of state, it is not possible to find an explicit bound on the couplings and the electric charge by applying the positivity constraints (3.12) and (3.13). However we can numerically investigate whether these physical constraints are satisfied whilst varying the cubic and quartic couplings for a given fixed charge. The corresponding pattern is given in figure 3, where we also check for positivity of the entropy (3.14) as well.

JHEP07(2019)012
The physical critical domain is the part of the parameter space for which the physical constraints discussed in section 3.2 are satisfied, with the property that a phase transition occurs. Looking at the left part of figure 3, for k = 1, this region has β 3 < 0 (or µ < 0) for all values of α 4 , with the exception of the axis β 3 = 0, where only for α 4 ≥ 0 are the associated phase transitions physical (the positive axis is green). The point β 3 = α 4 = 0 is also green, recovering the result that in the limit of vanishing cubic and quartic couplings, a charged AdS black hole still has physical critical points [52]. On the vertical axis α 4 < 0, we have γ 2 < 0. The entire physical region has only one critical point, whereas in the unphysical region it is possible to have either one, two, or three critical points depending on the given values of (β 3 , α 4 ). As we make the fixed value of e smaller (see for example the top right diagram in figure 3), we find that there are at most two possible critical points but only one of them is physical. Summarizing, in the presence of charge, critical points exist for all (β 3 < 0, α 4 ) (see figure 3). The center of the parameter space is the k = 1 Reissner-Nordstrom-AdS solution for which all constraints are satisfied in any dimension.

JHEP07(2019)012
Whenever critical points exist we find that the critical exponents 5 are α = 0, β = 1 2 , γ = 1, δ = 3, (4.15) which are the standard values from mean field theory, even when both numerator and denominator vanish in (4.13). These are typically obtained by considering the equation of state near the critical point [89], writing v = v c (φ + 1) , T = T c (τ + 1), (4.16) and expanding in powers of (φ, τ ). Since here we do not have a closed form for the critical quantities, we insert numerical values for parameters into equation of state to obtain critical values for T c and v c and we obtain yielding (4.15). We explicitly illustrate the occurrence of the phase transition by drawing a P − v graph in figure 4 for parameters for which there is a physical critical point by setting e = 1, β 3 = −4e 4 and α 4 = 5e 6 . We see clear Van der Waals behaviour, with two distinct phases for T < T c that coalesce at T = T c and become indistinguishable for T > T c . For sufficiently low temperature, the curve tends to negative pressures; however only positive values of the pressure are physical. The coexistence line is plotted in the right half of figure 4, illustrating the critical point at the end of a line of first-order phase transitions between large and small black holes.
The Gibbs free energy as a function of temperature is shown in figure 5, exhibiting the typical swallowtail characteristic of Van der Waals behaviour. It is notable that for quite small values of the pressure we still observe a swallowtail shape whose size grows rapidly. 5 The critical exponents quantify how physical quantities behaves in the vicinity of a critical point [52].
For t = T /Tc − 1, the exponent α characterizes the behaviour of the specific heat, while keeping volume constant The exponent β denotes a difference between the volume of a large black hole V l and the volume of a small black hole Vs on the isotherm process The behaviour of the isothermal compressibility κT is given by exponent γ The exponent δ characterizes the following difference on the critical isotherm T = Tc  Computing the specific heat we find that two stable branches of black holes exist, with the physical one at the global minimum of G.
For the allowed regions of parameter space in figure 3 for k = −1, the phase diagrams are qualitatively the same as in figures 4 and 5.

Critical behaviour in five dimensions
In this section we consider five dimensional solutions. The equation of state becomes 19) and the critical temperature is The critical volume satisfies the following equation and, as before, finding an explicit closed form for both the critical temperature and volume is not feasible. However for vanishing cubic coupling, there is a considerable simplification; solving (4.2) with β 3 = 0 for the corresponding critical temperature (4.20) yields we find where v c satisfies and we must have k = 1 so that t c > 0. This equation is a cubic polynomial in v 2 c and can be solved exactly. At most there are two real solutions for any given choice of parameters that are physically acceptable. 6 Figure 6 plots the number of critical points as a function of (β 3 , α 4 ) with fixed charge. For k = 1, unlike in 4 dimensions, we see that only if both couplings are non-zero we get two physical critical points in a certain region of parameter space, shown in dark green. The occurrence of two physical critical points for spherical black hole in five dimensions has to our knowledge not been seen previously. On the axes β 3 < 0, α 4 = 0 and β 3 = 0 (i.e., on the vertical axis) only for positive values of α 4 (or λ < 0) greater than some specific lower bound for the quartic coupling are the critical points physical, whereas for β 3 > 0, α 4 = 0 and β 3 = 0, α 4 < 0 they are unphysical, having γ 2 < 0. Physical critical points exist for JHEP07(2019)012 most of the region β 3 < 0, except for small values of |β 3 | and large enough large values of |α 4 | where γ 2 < 0. The case α 4 = 0 is discussed in [84].
Summarizing we have observed for the first time the occurrence of two physical critical points for spherical (k = 1) black holes and one critical point for hyperbolic (k = −1) black holes. We must have both couplings nonzero for this to take place.
Again we see that even if e = 0 there are physical critical points. For k = 1 we get regions with either one or two physical critical points, and only one for k = −1. In the latter case, there are some regions having negative mass (white region) and both couplings must be non-zero in order to get physical critical points. However in regions of parameter space having two physical critical points, there is new behaviour. We illustrate this in figure 7, which shows that there are two first order phase transitions. The transition at T = T c is standard VdW behaviour. But the second phase transition at T = T c is that of 'reverse VdW' behaviour: it is a transition from one phase for T < T c to two distinct small/large black hole phases for T > T c . Note that for a sufficiently large temperature the pressure becomes negative and the asymptotic structure of the spacetime is no longer AdS. Consequently there is an upper bound on the temperature of AdS black holes. We illustrate this in figure 7. Note that curves having P < 0 over a finite range of T can be given physical meaning via the equal-area law [92]. Referring to figure 7, we see that although P < 0 corresponds to a different asymptotic structure from that of AdS, the equal area law implies that P never actually attains these negative values, but rather remains constant and positive as the phase transition takes place. The coexistence line of the two distinct phase small/large has a critical point at a minimal value of T in contrast to that of a standard VdW phase transition. This phenomenon has been previously observed for black branes [34] and in cubic gravity [84]. In each plot, red lines depict the parts of the curves for which the specific heat is negative, blue lines indicate negative entropy, and the purple line indicates that both specific heat and entropy are negative. Top left: the low temperature region (for which there is a standard VdW phase transition), with P = 1.2P c (dotted, blue curve), P = P c (dotted, black curve) and P = 0.6P c , 0.2P c (solid, black and red curve) each plotted. The remaining graphs pertain to the high temperature region (for which there is a reverse VdW phase transition). Top right: P = 0.5P c , Bottom left: P = P c Bottom right: P = 0.99P c .

JHEP07(2019)012
The existence of a standard VdW phase transition followed by a reverse VdW transition at higher temperature is shown in figure 7. One might anticipate that with an appropriate choice of parameters that these two critical points (illustrated in the top right diagram of figure 7) could merge, yielding an isolated critical point (see the discussion for the six dimensional case in the next section). We find that this only occurs for either negative entropy and/or negative mass; these unphysical conditions persist in a neighbourhood of the merged critical point.
In figure 8, we illustrate the behaviour of the Gibbs free energy in both the lowtemperature region containing a standard VdW transition and in the high-temperature JHEP07(2019)012 region containing a reverse VdW transition. In the former case, we obtain the standard swallowtail behaviour for P < P c . In the latter case, for P sufficiently smaller than P c there are no phase transitions, as the upper right part of figure 8 indicates. As P → P c we obtain swallowtail behaviour, shown in the lower right part of figure 8.
For P = P c , in addition to negative specific heat we get regions with negative entropy (shown by the blue curve) and with both negative specific heat and negative entropy (shown by the purple curve); these are all unstable. For larger temperature the black hole solutions become stable (black line). Figure 8 shows that in addition to a reverse VdW transition, at high temperatures there are three black hole solutions with positive specific heat, with one having a minimal Gibbs free energy.
For the k = −1 hyperbolic black hole in d = 5 we get only a reverse VdW transition, in which higher temperatures have two distinct phases, as depicted in figure 9. The Gibbs free energy in figure 10 exhibits swallowtail behaviour for pressures a bit less than the critical pressure. However for larger values of pressure one of the branches corresponds to solutions with negative mass (depicted by the green curve) and the phase transition is no longer physical. For even lower pressure (bottom left diagram in figure 10) we see that unstable black holes with negative specific heat have lower free energy than those with positive specific heat. It is reasonable to expect that there is now a zeroth order phase transition between the two black curves in this figure. For even smaller P , the unphysical parts of the branches shrink and (apart from a small red region) become stable. As in the k = 1 case, there are again three black hole solutions with positive specific heat, with one having a minimal Gibbs free energy.

Critical behaviour in six dimensions
Turning now to six dimensions, we note from (2.6) that for the linearized field equations the contribution of the cubic term drops out 7 [24]. This yields some simplification, but the analysis is still somewhat complicated. The equation of state takes the following form

24)
7 In eight dimensions the quartic term drops out.

JHEP07(2019)012
with the explicit expression T c = 3 (36π 7 β 3 3 k + 24α 2 4 π 7 k)v c 18 + 180α 4 π 6 v c 16 β 2 3 k 2 + (162π 5 β 3 4 k + 276α 2 4 π 5 β 3 k)v c 14 +(256α 2 4 π 8 e 2 + 1476α 4 π 4 β 3 3 k 2 + 960π 8 β 3 3 e 2 − 144α 4 3 k 2 π 4 )v c 12 + (2880α 4 π 7 β 3 2 ke 2 for the critical temperature and being the equation yielding the critical volume. Setting the cubic coupling to zero, the associated critical temperature is 27) and the critical volume obeys the relation There is a singularity in the critical temperature at physical critical points; if α 4 < 0 we get γ 2 < 0 in the first case, and the critical temperature or critical volume becomes negative in the second case. We plot in figure 11 the possible critical points in the (β 3 , α 4 ) plane. Physical critical points appear only for β 3 < 0, and there can be as many as three for certain ranges of (β 3 , α 4 ) if k = 1 but only one for the hyperbolic (k = −1) case. Other possible critical points have one or both of γ 2 < 0 and S < 0. For vanishing charge we also have a region of at most two critical points for both k = 1 and just one critical point for k = −1; for vanishing α 4 there are still two physical critical points if k = +1, studied in detail in [84].
Clearly the maximal number of critical points for a given value of (β 3 , α 4 ) depends on both horizon geometry and dimension. We do not need to have both couplings non-zero JHEP07(2019)012 to obtain physical critical points. However in the absence of charge or when k = −1 (both with or without charge) both couplings must be nonzero.
The appearance of three physical critical points in the d = 6, k = +1 case stands in contrast to previous studies. This remarkable feature has not been observed before, and its occurrence is related to the conjunction of electric charge, cubic, and quartic couplings all being nonzero. Figure 12 shows that there is a reverse VdW transition in between two standard ones, one at cold temperatures T < T c 1 and the other at high temperatures T > T c 3 , with the critical temperature T = T c 2 of the reverse transition T c 1 < T c 2 < T c 3 .
The curves in figure 12 correspond to phase transitions that obey Maxwell's equal area law [92] with the actual pressure remaining positive during the phase transition (despite the curve indicating P becomes negative over a finite range of T ) as noted earlier. At sufficiently low temperatures the P − v curves cross the horizontal axis, yielding unphysical behaviour since the asymptotic structure is no longer AdS.
Critical points and coexistence lines are also displayed in figure 12. Numerical analysis confirms that for typical values of the parameters, each of the three critical points are characterized by mean field theory critical exponents, a hallmark of the end point of a first order phase transition. Note that the two critical points at high temperature are joined by a coexistence line.
For certain choices of the couplings and charge these two disjoint lines merge into each other, and an isolated critical point appears at the merge point. This new critical point is characterized by critical exponents that differ from the mean field theory ones. This phenomena was first observed in Lovelock and quasi-topological gravity [48,62,66], where isolated critical points were found to occur for hyperbolic horizons and massless AdS black holes. A thermodynamic singularity, at which pressure remains constant for any temperature and isotherms cross and reverse, was also found to occur at this point. However for Lovelock and quasi-topological black holes accompanied with conformal scalar hair, isolated critical points have been observed in five and higher dimensions for massive black holes without coinciding with a thermodynamic singularity [50,51]. More recently, in GQG cubic gravity only, isolated critical points have been found in six dimensions with no scalar hair for spherical black holes [84]. In all of these cases there are initially only two critical points that converge to one isolated critical point.
Here we see for the first time three critical points, two of which merge to form an isolated critical point. Furthermore these isolated critical points do not correspond to any thermodynamic singularity, i.e., the P − T curve does not have a zero slope. From the relation (4.17) we find that the coefficient B vanishes and the associated critical exponents arẽ which differ from the standard critical exponents appearing in (4.15) but are in accord with previous studies [48,62,66]. On either side of the isolated critical point the black holes satisfy all physical constraints, as the bottom right diagram in figure 12 indicates. The behaviour of free energy with respect to temperature for spherical black holes is illustrated in figure 13. At low temperatures, for P = P c the solution is stable and for P < P c a standard VdW transition occurs. For the two other critical points P c 2 , P c 3 there JHEP07(2019)012 is a region of the curve that has negative specific heat at low temperatures (red lines). For pressures P c 2 < P < P c 3 there are reverse and standard transitions that are shown by the rather sharp swallowtails depicted in the bottom graphs of figure 13 from left to right respectively. In contrast to the lower temperature (P < P c ) swallowtails depicted in the upper diagrams in figure 13, these swallowtails have positive specific heat everywhere apart from a few short segments. For regions of parameter space having two physical critical points, as depicted in figure 14 there is a first order standard VdW phase transition between small and large black holes at low temperatures and then a reverse VdW transition at higher temperatures. The right graph shows that for fixed charge and α 4 , and varying β 3 we obtain the appropriate value for the cubic coupling β 3 such that two coexistence lines meet, yielding an isolated critical point, with non-standard critical exponents given in (4.30). Again we see that there exist a range of temperatures on either side of the isolated critical point for which the black holes satisfy all physical requirements. Finally, in regions of parameter space with one physical critical point we get a standard VdW phase transition for k = 1. However if k = −1 then we find an reverse VdW transition that is similar what we described for the d = 5 hyperbolic black hole.

Critical behaviour in more than six dimensions
Increasing the value of d further, we find for seven dimensions that the qualitative features remain similar to d = 6 spacetime. Quantitatively, however, the parameter regions with only a single critical point get larger, whereas regions with two or three physical critical points get smaller. No further features emerge and so we shall not consider this case further.
However for d = 8 we find that we get up to three critical points. The structure of the associated phase transitions is similar to that we presented for d = 6 in the previous subsection, so again we shall not consider this case further.
Finally, we compute the ratio of critical quantities for any value of d. In general this must be done numerically, but we can obtain an analytic expression for small values of the couplings for which the physical constraints hold. Here we shall set the cubic coupling to zero 8 and compare results with the critical behaviour in Einstein gravity for spherical JHEP07(2019)012 black holes. To leading order the critical quantities read is the critical volume in Einstein gravity.
We see a dimension-dependent deviation from the critical values in Einstein gravity. The critical temperature and pressure increase and the critical volume decreases relative to Einstein gravity, except in four and six dimensions, where critical pressure and temperature are smaller and the critical volume is larger. In four dimensions these corrections are negligible.

JHEP07(2019)012
We finally obtain and we see that in the limit of vanishing quartic coupling, these results reduce those of charged k = 1 black holes in Einstein gravity [89]. The effect of the quartic curvature term on the Van der Waals ratio is to increase it in any dimension above the Einsteinian value. In four and five dimensions these corrections are negligible for small values of coupling and large enough charge. However in higher dimensions we see considerable deviation and the van der Waals ratio no longer is a "universal" quantity as it is in Einstein gravity: it depends on parameters of the theory as well as the dimension under considerations.

Thermodynamics in grand canonical ensemble
We now consider the grand canonical ensemble, in which we have fixed potential instead of fixed charge. From the viewpoint of AdS/CFT holography, fixed potential on the gravity side is related to fixed chemical potential on the CFT side. We first consider d = 4 and then discuss properties for generic dimensions, employing the approach in ref. [93]. We shall consider only the quartic term in the action; the cubic case was studied in [84]. We expect that when both couplings are nonzero a pattern similar to that of the fixed charge case for the number of critical points will emerge, though we shall not perform that analysis here.

Four dimensions
For fixed potential one needs to find expressions for the mass and the temperature solving again the equations of motion for the metric function near the horizon, since this choice of ensemble alters how the equations depend on the horizon radius. The first two leading JHEP07(2019)012 order terms in the expansion result in the following formulas, parameterized by the quartic coupling and r + : with the second equation yielding the equation of state where we defined v = 2r + . The explicit form of the critical quantities from the equation of state (5.2) are and the formula for the corresponding critical pressure is given inserting the above relations for T c and v c into equation (5.2). Note that we have the constraint so that critical quantities remain real. Solving equations (5.1) to determine the explicit form for M and T , and choose solutions that in the limit λ → 0 approach the Einstein branch, we find that only for k = +1 (spherical) black holes are there physical critical points.
Admitting a positive mass, while imposing the condition (2.6) with zero cubic coupling, we use the formula (4.8) to obtain an upper bound 0 < P ≤

JHEP07(2019)012
The black hole entropy at the critical point becomes and r + c = vc 2 with v c introduced in (5.3) and Z given in (5.4). Explicit numerical computation for k = +1 and λ < 0 indicates that the critical entropy is always positive; this will not hold for other values of these parameters.
In order to search for the phase structure of the black hole solutions, we obtain the formula for the free energy in the grand canonical ensemble which we plot in figure 15. We again observe standard Van der Waals behaviour. However we also see that the free energy is a decreasing function of the temperature for small T , vanishing at T → 0, in contrast to the fixed charge case; however a large branch of the curve has negative entropy. The P − v diagram in figure 15 illustrates that a phase transition happens for temperatures a bit less than the critical temperature. However for enough low temperatures, the pressure becomes negative, and we do not consider this unphysical case. Another way to observe the phase transition is via the coexistence line in the P − T plane, shown in the top right of figure 15. For low enough temperatures the entropy becomes negative, denoted by the blue solid line.
Instead of finding the equation of state for fixed potential, we perform the analysis with fixed pressure, commonly used in holography. The phase graph of temperature versus potential in figure 16 shows it again describes a first order phase transition between small and large black holes, and the coexistence line terminates at the critical point.

Higher dimensions
Here, we look into the thermodynamic properties of black hole solutions in the fixed ensemble in generic dimension. The equation of state in d dimensions is To compare results between the two ensembles with the cubic coupling set to zero, we note that in the fixed potential ensemble for spherical black holes in quartic GQG we get physical critical points for α 4 > 0 in four and five dimensions. However in six dimensions no physical critical points exist. This is in contrast with the fixed charge ensemble, where (α 4 > 0) we get single physical critical points in d = 4, 6 and two physical critical points in d = 5 (as well as single physical critical points).
In addition, for fixed potential and k = −1 hyperbolic quartic black holes, while there are potential critical points for d = 4, 6, these all have γ 2 < 0 (since α 4 < 0) and so there are no physical critical points for these dimensions (and likewise none for d = 5). This situation is the same as for the fixed charge ensemble, confirming that both cubic and quartic coupling must exist to obtain physical critical points.

Holographic hydrodynamics
One of applications of the AdS/CFT correspondence is the computation of the ratio of shear viscosity to entropy density η/s. We investigate in this section this ratio for the quartic theory.
It is well-known that for the field theories dual to Einstein gravity, the shear viscosity to entropy density ratio is lower-bounded by i.e., η/s = 1/(4π), and it has been suggested that this lower bound is universal, holding for any matter [94]. This conjecture thus states that η/s ≥ 1/(4π), and is the so-called KSS bound.
However higher derivative order contributions can cause the violations of this bound [95]. We therefore evaluate η/s for field theories dual to the quartic generalized quasi-topological theory for general d to see whether the KSS bound is satisfied.
Here, we employ the planar black hole solutions described by the following metric We define a new coordinate z = 1 − r 2 + /r 2 , to compactify the region outside the horizon. This gives where g(z) vanishes at z = 0, and g(1) = f ∞ . Expanding g(z) near the horizon 0 z 3 + · · · , (6.3) we solve the field equations to determine g (i) 0 for i = 2. As we discussed previously, the second derivative of the metric function near the horizon, g 0 , is not determined by the field JHEP07(2019)012 equations. However, its value can be chosen such that the numerical solution approaches its associated asymptotic solution. It is easy to check that by the coordinate transformation, the parameters g (i) 0 are written in terms of the parameters a i appearing in the near horizon expansion (2.28). We find where the explicit form of a 3 from section 2 for the planar black holes is Using methods described in [96], we perform a shift on the metric (6.2) with perturbation parameter . Computing the Lagrangian for the perturbed metric, and performing a series expansion for small , we obtain whereλ is the coupling in the action and is related to the coupling appearing in the equation of motion, using (2.3). The 'time formula' gives the shear viscosity Res z=0 √ −gL ω 2 2 , (6.8) whose explicit form for the case at hand is whereλ appears in the action (2.1) and is related to λ via (2.3).
Taylor expanding η/s about λ = 0 we obtain whereȧ 2 (0) denotes derivative of a 2 with respect to λ and then setting λ = 0. To compute the value of η/s numerically, one needs to determine the parameter a 2 for the choice of the other parameters, where the value of temperature is known from expansion near the horizon. In the above computation we replaced the parameter a 3 in terms of a 2 using the second order expansion near r = r + of eq. (2.10).
The Taylor expansion yields anȧ 2 (0) term that is left undetermined at this level. Only under the condition that the expression at first order in λ in the above series expansion for higher dimensions with some positive coupling (but in d = 4 for certain negative coupling) can we conclude that the KSS bound η/s ≥ 1/(4π) holds in all dimensions in the quartic generalized quasi-topological theories at small coupling.
The generalized quasi-topological term can cause the entropy density of black branes to change sign [34]. Hence there is a certain T p for which the ratio η/s exhibits a pole. Using the first order near-horizon expansion from (2.10), the corresponding quartic coupling is , (6.12) where the subscript "p" stands for "pole". It is interesting to note that in four dimensions λ p is equal the value λ c = −81 6 /1024 that we encountered in the critical limit of the theory in section 4.1; however, this coincidence does not hold in higher dimensions. However in any dimension it leads to the Einstein branch of the theory. The leading order term in (λ − λ p ) that describes the behaviour of the entropy near zero is

JHEP07(2019)012
Furthermore the corresponding expansion for the shear viscosity reads as 14) It is interesting to note that the entropy density vanishes linearly as λ → λ p . By contrast, the shear viscosity does not vanish in this limit as it contains a term consisting of some powers of a 2 evaluated at the pole. Explicit numerical evaluation shows that this term does not vanish as λ → λ p . Consequently the ratio of shear viscosity to entropy density has a pole at λ = λ p . Figure 17 shows that there is a smooth curve connecting η/s = 1/(4π) (for λ = 0) and η/s = ∞ (for λ = λ p ).
Since including quartic quasi-topological or Lovelock terms into the action does not alter the black hole entropy from its Einstein gravity value [34], only quartic GQG contributions are responsible for the occurrence of this pole. A previous study of cubic GQG found that similar behaviour for η/s took place [84].
To determine how the ratio η/s can be recast in terms of λ, we use the Padé approximant method to evaluate the parameter a 2 . For our consideration involving quartic term, the computations become cumbersome very rapidly while going to higher orders, so we only present only a few corresponding curves in four, five and six dimensions in figure 17. From this figure it is obvious that η/s starts from 1/(4π) for λ = 0 and then it grows until it diverges as expected at λ = λ p . For four dimensions this figure illustrates that the KSS bound for η/s holds.
From the relations (2.30) and (2.31) the explicit formulae for the temperature and mass are where we see that T is independent of λ in five dimensions. In figure 17 we plot η/s for the values of λ such that physical constraints are satisfied.
To determine this, one needs f ∞ > 0 for an asymptotic AdS spacetime and P (f ∞ ) > 0 to satisfy the no-ghost criterion. The behaviour of f ∞ in four dimensions is given in right diagram in figure 2; similar behaviour is observed in d = 5, but with slightly different contributions from the different solution branches. In four and five dimensions we find that λ must be negative so that γ 2 > 0. The expression for the mass in (6.15) sets a lower bound on the coupling. In conjunction with the condition γ 2 > 0 in (2.23) required for a well-defined asymptotic region, we must have 0 > λ > λ p = − 6 /16 for any choice of AdS length and black hole horizon, where the lower bound λ p corresponds to zero mass and vanishing entropy; note that |λ p | is smaller than |λ c | given in (4.9). We note from figure 17 that the KSS bound is violated for λ/λ p > 0 or in other words for all physically acceptable values of λ since λ p < 0.
In six dimensions solving equations (2.30) and (2.31) yields four different solutions for the mass and temperature. However the only branch that satisfies the criteria for a physical solution is Here, the mass is positive for 0 < λ < λ p and vanishes (along with the entropy) at λ = λ p ; from eq. (2.24) we find

JHEP07(2019)012
which is positive both for λ > 0 and for sufficiently negative λ. On the other hand, even though the overall structure of vacuum solutions looks like what is exhibited in the right graph in figure 2, in this case the upper branch with negative λ is associated with the existence of ghosts -it is therefore excluded from further consideration. The lower negative λ branch yields negative γ 2 . Hence only positive values of λ yield physical solutions as well as γ 2 > 0. However there is an upper bound for the positive coupling that is enforced by positive mass. From figure 17 starting from the same value of the ratio at zero coupling, for small absolute values of λ, in four dimensions the ratio is initially larger than for d = 6 but then declines to smaller values. Both curves blow up as the pole is approached. In these cases the KSS bound holds, and mass remains positive as 0 < |λ| < |λ p |.
We anticipate similar behaviour in seven and higher dimensions, namely that positive coupling is required for correct asymptotic behaviour and physical mass, and also satisfies the KSS bound.
It is known that there is an upper bound on the coupling restricting the existence of acceptable CFT duals, so the whole range of λ ∈ (0, |λ p |) cannot possess a holographic interpretation [97,98]. We postpone how to address this issue to future investigations.
If both cubic and quartic couplings are nonzero, there is still a point in parameter space for which the entropy vanishes. This is a singular point. It cannot be removed since a 2 , which appears at zeroth order in the η expansion (similar to eq. (6.14)), is computed for the values of the couplings at the pole. Corresponding results for cubic gravity are discussed in [84].

Discussion
We have investigated charged static spherically symmetric AdS black holes for both spherical (k = 1) and hyperbolic (k = −1) geometries in generalized quasi-topological gravity (GQG). These theories are of notable interest since this class of solutions has a single metric function, analogous to Lovelock and quasi-topological gravity at the same order. We have considered both cubic and quartic GQG to see how these additional terms modify the results for Einstein gravity in four, five, and six dimensions.
Although the metric function cannot be obtained analytically, it is feasible to find both the asymptotic and near horizon behaviour of the metric perturbatively. We then apply the shooting method to verify that these solutions match in the intermediate region. The near horizon expansion characterizes the mass and temperature of black holes, and therefore, the thermodynamics of the black hole can be completely understood, despite the lack of an exact solution. Furthermore our numerical considerations demonstrate that for either fixed cubic coupling or quartic, increasing the electric charge correspondingly decreases the horizon radius. On the other hand at fixed charge, enlarging the coupling has the effect of inceasing the horizon radius. Investigating the solutions near the origin r = 0 in four dimensions we find that the curvature scalar singularity is softened.
Taking the cosmological constant and cubic and quartic couplings as thermodynamic variables, we examined the thermodynamic properties in the given spherically symmetric JHEP07(2019)012 configuration in detail, including verification of the extended first law and Smarr relation. However not all solutions obey standard physical requirements (positive mass, positive entropy, AdS asymptotics, and the no-ghost condition), so we constructed the physical constraints between the couplings and the charge that gives the domain for parameters to yield physical critical points.
In some regions of parameter space black hole entropy can be negative. The sign of the entropy depends on the spacetime dimension and is a function of temperature in terms of horizon radius. For a fixed charge ensemble or for neutral black holes, there are situations (for small black holes) in which the entropy becomes negative as r + → 0. One can add the absolute value of this amount to the entropy to ensure that the vacuum state has zero entropy. Under these circumstances the associated free energies are shifted by the same amount, so the thermodynamic properties are not affected.
Another way to shift entropy to a positive value is by adding an explicit Gauss-Bonnet contribution to the action [32,99] or by including the volume form of the induced metric on the horizon in the Lagrangian [100] where, similar to adding the absolute value of the negative entropy as r + → 0 as mentioned above, one must ensure S → 0 as M → 0.
More generally, in even dimensional spacetimes adding Euler densities to the action will shift the entropy by an arbitrary constant without changing the solution of the field equations (see appendix A of [84] for the case of Gauss-Bonnet gravity in four dimensions). This arbitrary constant must be chosen such that the entropy vanishes for AdS spacetime. The more general problem of how to deal with negative entropy in any dimension we leave to future work.
Working in both the fixed charge ensemble and fixed chemical potential ensemble, we classified the phase structure and critical points for these black holes. In the fixed charge ensemble, and in four dimensions we found out that even for zero charge, there are still physical critical points for k = 1 when at least one of couplings is non-zero (in contrast to Einstein gravity) and also for k = −1 when both couplings are nonzero. A first order VdW phase transition between small and large black holes was seen.
For the first time we have observed critical points for a neutral hyperbolic black hole in any dimension provided both the cubic and quartic couplings are nonvanishing. This emphasizes the importance of the non-linearity of GQG in inducing new phase behaviour.
In five dimensions we observed the occurrence of two physical critical points when both cubic and quartic couplings are nonzero, even for vanishing charge. These respectively correspond to the end point of a standard VdW transition and the starting point of a reverse VdW transition. However for the reverse VdW transition the pressures are smaller than the second critical pressure, and so one cannot choose parameters such that these two end points merge to obtain an isolated critical point -the pressure in the second line of firstorder phase transitions does not increase with temperature unlike the first coexistence line.
For hyperbolic black holes in five and six dimensions there are regions in parameter space that yield negative mass. In the d = 6, k = +1 case, we noted the existence of three and two critical points under the respective condition that three and two of the parameters (charge and two couplings) are non-zero. Since the coexistence lines are both increasing JHEP07(2019)012 functions of pressure with respect to temperature, it is possible to find parameter choices for which the critical points merge into an isolated critical point.
We obtained the critical temperature and volume in terms of the electric charge and couplings in various dimensions up to first order in the quartic coupling λ. For spherical black holes in quartic gravity the universal relationship in Einstein gravity given by the first term in (4.32) receives dimension-dependent corrections -it is no longer universal. These correction are negligibly small in four and five dimensions; however in higher dimensions they cause the critical ratio (4.32) to increase.
We also analyzed the existence of physical critical points in the grand canonical ensemble in the quartic theory. Obtaining the relevant thermodynamic quantities, we confirmed the presence of a first order phase transition in four dimensions that is absent in the corresponding situation for Schwarzschild black holes in both cases, for which either the chemical potential or the pressure are fixed. We also investigated phase transitions in higher dimensions.
Our study of black hole thermodynamics will, of course, be modified by the inclusion of Lovelock terms in d > 4 and by standard quasitopological terms [16,22,101]. While it is conceivable that such terms in combination with those in the quartic theory will yield new phase behaviour, their inclusion would considerably broaden the parameter space and so we shall leave this investigation for future work.
Our goal was to isolate the thermodynamic behaviour of black holes in the quartic theory. The reason for omitting the Lovelock terms (which in principle could be present in more than 4 dimensions) was primarily for simplicity. Their inclusion broadens the parameter space, and would make our paper notably longer than it already is. A study of how the Lovelock terms affect our results is possible, of course, but we prefer to leave this for future work.
Finally, in the context of the AdS/CFT, we computed the ratio of shear viscosity to entropy density η/s for field theories dual to the quartic generalized quasi-topological theory in all dimensions and concluded that in four dimensions the KSS bound held for choices of the quartic coupling yielding positive mass and temperature. In striking constrast, we found that this behaviour did not hold in five dimensions -the range of coupling required for a physical solution entailed a violation of the KSS bound. However, the bound remains valid in four and six dimensions, and we anticipate it is satisfied in higher dimensions. It would be interesting to see if including yet higher order terms in the curvature can yield black holes satisfying both the KSS bound and all physically reasonable requirements in five dimensions as well as other dimensions. Further constraints could be imposed by configuring other conditions in the corresponding CFT, such as causality constraints and positivity of energy flux; the implications of this require further investigation.

JHEP07(2019)012
The near horizon expansion for the metric function is f (r) = 4πT (r − r + ) + a 2 (r − r + ) 2 + ∞ i=3 a i (r − r + ) 3 . (B.1) As discussed earlier the field equations at each order determine the parameters in the expansion in terms of a 2 , but a 2 = f (r + )/2 itself remains undetermined. One approach for fixing this free parameter is to use the shooting method such that the numerical solution for the metric function reaches the known asymptotic solution at large r. Another method [9] entails considering a 2 = g(λ).
Continuing to higher orders, we find that g (n) (0) is specified at any arbitrary order provided a n+3 does not have a singularity as λ → 0, and a n+2 receives the expected Einstein value at the same limitation. Because the coefficients of the derivatives increase rapidly, the Taylor series does not have a non-zero radius of convergence and the analytic expression for a 2 is not valid. However inserting the g (n) (0) terms (which implicitly contain derivatives of the temperature) into a Padé approximant yields a satisfactory result. The coefficients of the Taylor series diverge due to the existence of a pole at positive λ. This problem comes from the fact that in our calculations we consider the temperature as a function of the coupling, although it is not a real analytic function.
Computation of the Padé approximants for higher orders is quiet tedious, although according figure 18, for small coupling one can make use of lower order Padé approximant to obtain accepted outcome. For example, [2|2] Padé approximants for some dimensions in terms of x = λ/λ p are, a d=4 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.