Extremal Rotating Black Holes in Einsteinian Cubic Gravity

We obtain new solutions of Einsteinian cubic gravity coupled to a Maxwell field that describe the near-horizon geometry of charged and rotating black holes. We show that the AdS$_2\times\mathbb{S}^2$ near-horizon geometry of Reissner-Nordstr\"om black holes receives no corrections, but deviations with respect to the extremal Kerr-Newman solution appear as we turn on the angular momentum. We construct the profile of these corrected geometries using both numerical methods and slowly-spinning expansions, but we also find additional solutions that do not reduce to AdS$_2\times\mathbb{S}^2$ geometries in any limit and that do not have a counterpart in Einstein gravity. Quite remarkably, we are able to obtain closed-form exact expressions for the area and Wald's entropy of all of these solutions, and using this result, we analyze the phase space of extremal back holes in this theory. To the best of our knowledge, this is the first time the entropy of a rotating black hole in higher-order gravity has been exactly computed.

We obtain new solutions of Einsteinian cubic gravity coupled to a Maxwell field that describe the near-horizon geometry of charged and rotating black holes. We show that the AdS2 ×S 2 near-horizon geometry of Reissner-Nordström black holes receives no corrections, but deviations with respect to the extremal Kerr-Newman solution appear as we turn on the angular momentum. We construct the profile of these corrected geometries using both numerical methods and slowly-spinning expansions, but we also find additional solutions that do not reduce to AdS2 ×S 2 geometries in any limit and that do not have a counterpart in Einstein gravity. Quite remarkably, we are able to obtain closed-form exact expressions for the area and Wald's entropy of all of these solutions, and using this result, we analyze the phase space of extremal back holes in this theory. To the best of our knowledge, this is the first time the entropy of a rotating black hole in higher-order gravity has been exactly computed.

I. INTRODUCTION
General Relativity describes accurately the dynamics of the gravitational field in the regime of relatively low curvature, but modifications of this theory are expected to appear at high energies. The fact that GR is incompatible with quantum mechanics [1][2][3] indicates that it should be regarded as an effective theory, presumably arising from an underlying theory of quantum gravity. Independently of what the UV-completion of GR turns out to be, it is broadly accepted that an effective low-energy description of that theory will contain the Einstein-Hilbert action plus an infinite tower of higher-derivative corrections -this is, in particular, a definite prediction of String Theory [4][5][6][7][8][9][10]. Such corrections modify the behaviour of the gravitational field when the distances involved are of the order of the length scale of new physics. Thus, they become extremely relevant in the very early universe or near black hole singularities, but also at the level of the horizon of small enough black holes. It is therefore an interesting task to determine the properties of the modified black hole solutions, with particular emphasis on the corrections to the thermodynamic quantities, such as entropy and temperature [11][12][13][14][15].
From the point of view of Effective Field Theory (EFT), one should treat the higher-derivative corrections as perturbations over the GR geometry. Obtaining the corrected solutions in this perturbative approach is usually an accessible task; however, perturbative solutions give us very little information. In fact, the perturbative corrections are only valid as long as they remain very small, and many potentially interesting phenomena, that would appear at a non-perturbative level, are lost. For this reason, it is interesting to find exact black hole solutions of higher-order gravity. * pabloantonio.cano@kuleuven.be † david.perenniguez@uam.es The problem of obtaining exact black hole solutions is, of course, more complicated. Let us consider first the case of spherically symmetric black holes. Until very recently, the only theories in which exact solutions modifying in a non-trivial way the Schwarzschild geometry had been constructed were Lovelock [15][16][17][18][19][20][21][22][23] and Quasi-topological gravities [24][25][26][27], both types of theories existing only in D > 4 dimensions. 1 The gap in D = 4 has been recently filled thanks to the construction of a new type of theories with very interesting properties. Known as Generalized Quasi-topological gravities (GQTGs) [32], these theories allow for simple spherically symmetric black hole solutions whose thermodynamic properties can be studied analytically [32][33][34]. Besides, GQTGs exist in all dimensions (including, in particular, D = 4) and at all orders in curvature [35], and very likely they provide a basis to construct the most general EFT for gravity [36]. Spherically symmetric solutions in these theories have been studied at all orders in curvature in D = 4 [37] and at cubic [32,38] and quartic order [34,39] in various dimensions, and this has allowed us to gain substantial information about spherically symmetric black holes in higher-order gravity. In particular, one of the most remarkable features of these theories is that black holes become stable below certain mass [37], hence avoiding the complete evaporation in a finite time and the final explosion of black holes. This is analogous to the behaviour of higher-dimensional Lovelock black holes found long ago in Ref. [11]. In this paper we will consider an extension of Einstein gravity containing the simplest nontrivial Generalized Quasi-topological density in D = 4, which is known as Einsteinian cubic gravity (ECG) [40]. This theory was the first member of the GQT class to be discovered and we review some of its properties as well as recent results in Sec. II.
Despite the success in the construction of spherically symmetric black holes, a remaining issue in the world of higher-order gravities is to find rotating black hole geometries. 2 In fact, exact rotating solutions have not even been found in Lovelock gravity, which is the simplest nontrivial extension of GR that one could consider. 3 Thus, the question about what a rotating black hole in higherderivative gravity is like has not been answered yet. However, this is a primordial question, since, after all, realistic black holes will in general possess angular momentum.
The equations of motion for an axisymmetric and stationary metric are far more complicated than those in the spherically symmetric case. Even though we expect some simplification of the equations taking place for GQTGs -because they do so in the static case -, obtaining a complete rotating black hole solution would necessarily require a laborious numeric computation. However, there are several limits in which the problem is simplified. On the one hand, one might consider slowly-rotating solutions and stay perturbative in the spin. This has been explored in the case of quadratic [44] and cubic [45] Lovelock gravity. The case for D = 4 ECG will be reported in a coming publication [46]. On the other hand, it is possible to study the opposite limit, namely, the case of extremal black holes. In this situation, the horizon is placed at an infinite distance and the near-horizon limit is well-defined, giving rise to a new solution of the gravitational equations. This near-horizon geometry has more symmetries than the global solution, and this enormously simplifies the problem of solving the field equations. We will show in this paper that the equations of motion of ECG reduce in this case to a single secondorder ODE. This equation has to be solved numerically, but most remarkably, we will see that it is possible to obtain the exact expressions for the area and entropy of these black holes without using any approximation. We are not aware that a similar analysis has been performed for other pure-metric higher-order gravities, but let us mention that Ref. [47] computed the (perturbative) corrections to the near-horizon geometry of extremal Kerr black holes in the case of Einstein-dilaton-Gauss-Bonnet [48][49][50] and dynamical Chern-Simons [51] gravities.
For generality purposes, we will add as well a Maxwell field into the game, which will allow us to study rotating and charged extremal black holes. This will prove to be useful, as AdS 2 × S 2 geometries -corresponding to nonrotating charged black holes -are always solutions of higher-order gravities. The rotating black holes can then be studied as a deformation of these geometries, which facilitates the analysis of the solutions. However, we will also show that there are new branches of solutions that do not reduce to AdS 2 ×S 2 geometries in any limit. These solutions do not exist in the Einstein gravity limit and, as we will see, they have somewhat exotic properties.
The paper is organized as follows. We start in Sec. II by introducing our theory, corresponding to ECG coupled to a Maxwell field. In Sec. III we write the metric and vector ansätze for a rotating near-horizon geometry possessing an SL(2, R) × U(1) isometry group, and we evaluate and partially solve the equations of motion. We reduce the field equations to a single second-order ODE for one variable. Then we discuss the boundary conditions that need to be imposed in order to obtain fully regular solutions. In Sec. IV we study in detail the solutions of the previous equation that are smooth deformations of AdS 2 × S 2 geometries. We construct solutions -both numerically and in the slowly-rotating approximation -which are labeled by the total charge Q and by a parameter x 0 which we argue is related to the spin a = J/M . More interestingly, we find that both the area and the Wald's entropy can be obtained exactly, and we study them as functions of Q and x 0 . In addition, the physically meaningful relation S(A, Q) is derived and we also study its profile. In Sec. V we analyze the full space of near-horizon geometries, showing that there exists an important degeneracy of solutions. We discuss the properties of the additional branches and comment on the structure of the diagram S(A, Q). Finally, we draw our conclusions in Sec. VI. We also include a number of appendices that support and extend some of the results in the main text.

II. EINSTEINIAN CUBIC GRAVITY
We consider the following theory which consists of the (cosmological) Einstein-Maxwell action -where F µν = 2∂ [µ A ν] -plus a cubic curvature correction P, the Einsteinian cubic gravity density [40] Also, µ is a dimensionless coupling while L is a length scale that determines the distance at which the gravitational interaction is modified. As stated earlier, P is the lowest-order non-trivial member of the GQT family of theories in D = 4. On a historical note, this theory was first identified by the special form of its linearized equations on maximally symmetric backgrounds, which turn out to be of second order in any dimension [40]. Afterwards, the simple form of spherically symmetric black hole solutions in this theory was noticed [52,53], and this triggered the construction of the GQT class of theories [32][33][34]. By now, many other aspects of ECG have been explored, including the characterization of observational deviations with respect to GR [54,55], holographic applications [56][57][58], inflationary cosmologies [59][60][61] 4 and other types of solutions [62][63][64].
Up to the six-derivative level, P represents the leading parity-preserving higher-derivative correction to the Einstein-Hilbert action [36]. However, when a Maxwell field is included, there are other terms that we could add at this order. Schematically, these would be of the form F 4 , RF 2 , F 6 , RF 4 , R 2 F 2 . Nevertheless, it is not our intention to study the most general correction to extremal Kerr-Newman geometries. Instead, we focus on the theory above because it will allow as to perform many explicit computations that are practically unaccessible for other higher-derivative theories.
The equations of motion of (1) read where the gravitational tensor E µν and the energymomentum tensor T µν are given by

III. NEAR-HORIZON GEOMETRIES
Near-horizon geometries of extremal rotating black holes possess an isometry group SL(2, R) × U(1), and a general ansatz for this type of metrics can be written as 4 In the cosmological context, the solutions appearing in Refs. [59][60][61] were constructed in a modified cubic theory that takes the form P − 8C, where P is the ECG term -see (2) -and C is a cubic piece that makes no contribution on spherically symmetric backgrounds [52]. Thus, that combination produces the same black hole solutions as the ones constructed for ECG. We have checked that C is irrelevant for our present setup too. [65] which depends on two functions f (x) and N (x) and on one constant n. In addition, we consider a vector field of the following form which depends on another function h(x). Then, we have to insert this ansatz in the equations of motion (3) and solve them. Due to the symmetries of the ansatz, one can check that the only independent components of the Einstein's equations are those corresponding to µν = xx and µν = ψψ -the rest are related to them by the Bianchi identities. Thus, we only need to solve those equations together with the Maxwell equation.
An important observation is that these equations allow for solutions that have N (x) = 1. The reason is that, when evaluated on N (x) = 1, the components of the gravitational tensor -which we show in Appendix Abecome proportional, namely and the same property holds for the Maxwell energymomentum tensor T µν . In general, higher-derivative gravities do not satisfy the condition (9), meaning that these theories do not allow for solutions with constant N (x). In turn, it is quite remarkable that this property holds for ECG. As we are going to see, this represents a drastic simplification of the equations of motion. Let us also note at this point that, besides the solutions with N (x) = 1, there can be other solutions. In fact, Einstein gravity allows for solutions with non-constant N (x), but these turn out to be singular, and only the solutions with N (x) = 1 represent the near-horizon geometry of extremal Kerr-Newman black holes. In the same way, ECG will presumably allow as well for this type of singular solutions when N (x) is non-constant. Thus, from now on we set N (x) = 1. Now, we can evaluate Maxwell's equation, which turns out to be independent of f (x): where a and b are two integration constants that are related to the electric and magnetic charges. Thus, at this point we have reduced the problem to solving one equation for f (x), namely E xx = T xx . However, before going into the resolution of this equation, let us massage a bit the solution in its current form. Let us note that the coordinate x is compact and it will range within two values x 0 > 0 and −x 0 . These values are determined by the vanishing of the function f (x) -which is assumed to be even -at those points: f (x 0 ) = f (−x 0 ) = 0. Also, let us introduce the quantity Then, observe that in order to avoid a conical singularity at x = ±x 0 -these points will correspond to the poles of the horizon -the coordinate ψ must have period 2π/ω. Using these results, we can already compute the electric and magnetic charges even if we do not know explicitly the function f (x). In Planck units, these charges read where the integration is performed on any surface of constant t and r. We note that these are the actual values of the charges that we would obtain in a global solution containing an asymptotic region. Let us finally exchange x and ψ in terms of two new coordinates (17) and let us introduce the function In this way, we rewrite our solution in the following form and by construction, g(y) satisfies Let us finally evaluate the remaining equation, which in the new coordinates is E yy = T yy . On the one hand, we have where Q 2 = q 2 + p 2 . On the other hand, E yy takes the form of a total derivative, namely where Hence, integrating both sides of the equation we obtain where N is an integration constant. Thus, we have reduced the equations of motion to a single ODE of second order for g(y).
Our task now is to solve the previous equation in or-der to obtain near-horizon geometries. So far, we have included a non-vanishing cosmological constant for generality, but for the sake of simplicity we set Λ = 0 from now on. The case of Λ = 0 is briefly discussed in Appendix D.

A. Einstein gravity
Let us first of all check that we recover the near-horizon geometry of extremal Kerr-Newman black holes when we set µ = 0. In that case, Eq. (25) is simply algebraic and we obtain the solution straightforwardly, . (26) We can see that the parameter N breaks the symmetry y ↔ −y of the solution that we assumed in identifying the charges q, p. More importantly, when N is present (and x 0 = 0), there is necessarily a conical singularity at one of the poles of the horizon (where g vanishes), because the slope of g will be different in each one. In fact, N is the NUT charge, and it is known that NUTcharged, rotating black holes present this type of conical singularities at the horizon [66]. In order to avoid these problems, we set N = 0. In that case, g(y) is even, and we have to impose the conditions (21), which are going to fix several relations between the parameters of the solution. We find and after simplifying we obtain We see that this is the near-horizon geometry of extremal Kerr-Newman black holes (NHEKN) [67], where x 0 is nothing but the angular momentum per mass x 0 = a. Likewise, n = M is the total mass and ω is the angular velocity of the horizon. In addition, we can compute the area, which reads For x 0 = 0 we recover AdS 2 × S 2 , which is the nearhorizon geometry of extremal Reissner-Nordstrom black holes.

B. Einsteinian cubic gravity
Let us now consider a non-vanishing µ. In analogy to the Einstein gravity case, we set the NUT charge to zero, N = 0, in order to avoid conical singularities. Now, once the corrections are included, the equation (25) becomes of second order and we need to impose appropriate boundary conditions in order to solve it. We warn that the constraints (21) are not really boundary conditions: they are restrictions to the parameters of the solution. Instead, the boundary conditions we will impose are the following: (1) the solution is even, and this is equivalent to asking g (0) = 0. (2) The solution is regular at y = ±1, i.e., it is analytic at those points. Therefore, according to (21), the solution should have a Taylor expansion near y = 1 of the form for some coefficients g k . When this expansion is inserted in (25) we can Taylor-expand the equation as well, obtaining the following series (31) Thus, all the coefficients C k must vanish and this gives us a series of equations for the parameters of the solution.
Remarkably, the first two equations C 0 and C 1 are independent of the g k , and instead they provide two relations between x 0 , n, ω and Q: We have seen that in Einstein gravity x 0 is identified with the angular momentum per mass, a, while in turn n is the mass and ω is the angular velocity. We cannot expect that the same identifications work for higher-curvature gravity, and, since we lack the asymptotic region, we cannot correctly identify these quantities. Nevertheless, since x 0 controls the degree of non-sphericity of the solution, we do expect that there will be a monotonous relation between this parameter and the angular momentum -we recall that this parameter enters in the metric as ds 2 = (x 2 0 y 2 +n 2 )ds 2 AdS2 +. . .. Hence, it seems reasonable to use x 0 and the charge Q to label our solutions. Then, the equations (32) provide us with the values of n(x 0 , Q) and ω(x 0 , Q). It is worth emphasizing that such equations are exact; we have implemented no approximation in our approach. Besides, this allows us to compute the area of these black holes even if we do not know g explicitly, since it is given by Then, once the parameters n and ω (or alternatively A) are determined, we can solve the rest of the equations C 2 = 0, C 3 = 0, etc. It tuns out that these equations fix all the coefficients g k>3 in (30) in terms of g 2 , which is the only free parameter. Thus, we find that there is only a one-parameter family of solutions that are regular at the pole y = 1, which means that regularity is in fact fixing one integration constant. Now, the remaining parameter g 2 is determined by the condition that g be an even function, which is equivalent to asking g (0) = 0. Thus, we have a well-defined boundary problem, which at most will possess a discrete number of solutions.
Let us summarize our findings so far. Our near-horizon geometries are labelled by two parameters which we can choose to be Q and x 0 . Imposing regularity of the solution at y = ±1 yields the equations (32), whose solutions give the possible values of n and ω. Finally, the differential equation (25) must be solved imposing the regularity condition (30) and g (0) = 0. As we will see later, the equations (32) have more than one solution for fixed Q and x 0 , which leads to an important degeneracy of near-horizon geometries that have the same Q and x 0 . However, it turns out that there is only one branch of solutions that are smoothly connected to an AdS 2 × S 2 geometry in the limit of x 0 → 0. In this section we will focus our attention on those solutions.
Let us first solve the equations (32) when x 0 << Q by assuming a series expansion of the form ω = n ω k x k 0 , n 2 = k c k x k 0 . We find the following solution where the higher-order terms can be easily computed as well and we show few of them in Appendix B. Now, let us also assume a series expansion of the metric function g(y), so that Plugging this expansion together with (34) in the equation (25) we find the equation satisfied by every component g k (y). The leading term g 0 -which is the only one that survives in the limit x 0 → 0 -satisfies the following equation We can see that a solution of this equation fulfilling the appropriate boundary conditions is given by Thus, in the limit x 0 → 0 the metric (19) becomes which corresponds to an AdS 2 × S 2 geometry. In fact, this is the near-horizon geometry of extremal Reissner-Nordstrom black holes, and, as we can see, it possesses no corrections. Thus, this is an exact solution of ECG for any value of µ. Let us then consider the effect of rotation by assuming a finite x 0 . Analyzing the equations for the following terms, g k (y), we see that they all allow for a solution which is a polynomial in y, and that this solution is the only one that satisfies the boundary conditions. For instance, up to quadratic order in x 0 we have and more terms are shown in the Appendix C. A few comments are in order. First, let us remark that this is a perturbative expansion in x 0 , but it is exact in µ. Second, we observe that if we put µ = 0 in the expression above we get g(y) , which coincides with the perturbative expansion of the NHEKN solution (28), and the same holds for the higherorder terms that we show in the appendix. Therefore, these solutions in principle approach the NHEKN one when µ → 0, or more precisely, when Q >> µ 1/4 L, i.e., when the size of the black hole is much larger than the length scale of the corrections. However, there is a sub-tlety: we observe that the O(x 2 0 ) term (and also all the higher-order ones) diverges for Q 4 = 9µL 4 . In general, we observe that all the terms of order greater or equal than 2n diverge for Q 4 = 3((n+1) 2 −1)µL 4 . This implies that, when Q crosses one of these values, the solution changes discontinuously, and near those critical values we seem to find no solution. Therefore, as we increase Q and x 0 , the solution will approach the NHEKN one, but it will make it in a non-continuous way. This is best understood by constructing the nonperturbative numerical solutions. We show several of them in Fig. 1, where we represent the function h(y) ≡  ω(1−y 2 ) , which allows for a more direct comparison between the different curves. We have checked that, when x 0 is small enough, the slowly-rotating expansion (39) gives a very good approximation to the numerical curves. Looking at Fig. 1 we observe that, indeed, the profile of the solution is quite different for distinct values of Q, but eventually it becomes similar to the NHEKN one for large black holes. In addition, in Fig. 2 we show the embedding of the black horizon in Euclidean space for some of these solutions.
One important drawback, though, is that we do not seem to find solutions when x 0 /Q is large. As we can see, Eq. (25) becomes singular at the points in which x 4 0 yg = 0, which implies that the coefficient of g vanishes. This only happens when the ratio x 0 /Q is large enough. For example, if we evaluate the previous expression for NHEKN geometries and we ask that it does not vanish at any point, we must impose x 0 /Q < 1/ √ 3. Now, if that quantity vanishes, the solution will typically become singular at that point, unless we fix a regularity boundary condition there. But in that case, we cannot impose the boundary conditions of regularity at y = ±1 and that g (0) = 0. Hence, we find that, even in the regime where the corrections are small, the equation (25) has no regular solutions correcting the NHEKN geometry for x 0 /Q large.
In addition, our numerical exploration indicates the existence of an important multiplicity of solutions even when the boundary conditions are fixed. This is, once we have solved (32) and found n(x 0 , Q), ω(x 0 , Q), the equation (25) seems to have different solutions that differ on the profile of g(y). This already happens in the x 0 → 0 equation (36), which possesses other solutions than (37) satisfying g(±1) = 0, g (±1) = ∓2/Q 2 . These do not need to be similar to the NHEKN geometry even when µ is small, and in general they will possess a different domain of existence from the solutions considered in the preceding paragraph. In any case, all of these solutions are characterized by the same set of parameters x 0 , Q, n, ω, so they share a number of common properties.
In order to simplify the discussion, in the next subsection we will remain agnostic about the existence or non-existence of solutions of Eq. (25). Providing some solution exist, we are going to see that the area and entropy can be obtained exactly without knowing the profile of g(y).

A. Area and entropy
As we have seen, it is possible to solve the equation (25) either perturbatively in x 0 or numerically. Nevertheless, there are some properties of these near-horizon geometries that we can compute exactly. One of them is the area, which is given by (33). Then, using Eqs. (32) one is able to obtain the area as a function of x 0 and Q. The relation A(x 0 , Q) for several values of Q is shown in Fig. 3. Near x 0 = 0, one can use the expansions (34) in order to obtain the approximation which is valid as long as x 0 << Q. Thus, for x 0 → 0, the area reduces to the corresponding value of extremal Reissner-Nordstrom black holes, but looking at Fig. 3 we see that an interesting behaviour takes place when we increase x 0 . If the charge is large enough, the corresponding curve differs slightly from the value in Einstein gravity for intermediate values of x 0 , but for large x 0 one recovers again the extremal Kerr-Newman result A ∼ 4π(Q 2 + 2x 2 0 ). On the other hand, if the charge is too small -the threshold value is -the curve does not approach the Einstein gravity result, and instead we see that A tends to a constant for x 0 → ∞. This represents an exotic solution that does not exist in Einstein gravity, and it satisfies where α is a constant determined from the equation On the other hand, near-horizon geometries allow us to compute the entropy of black holes, even if we do not know the behavior in the asymptotic region, thanks to Wald's entropy formula [12][13][14], which reads 5 In this expression, the integral is taken over the horizon H, h is the determinant of the induced metric on H and µν is the binormal of the horizon, normalized as µν µν = −2. Applying Wald's formula (45) to our theory (1), we get where P µναβ is the tensor defined in (6). The horizon of the metric (19) is placed at r = 0, but the integration can be equivalently performed on any slice of constant t and r. The non-vanishing components of the binormal read tr = − rt = (y 2 x 2 0 +n 2 ), so that P µναβ µν ρσ = 4(y 2 x 2 0 + n 2 ) 2 P trtr . Remarkably, we find that this quantity takes 5 For Lagrangians containing covariant derivatives of the Riemann tensor, the partial derivative of the Lagrangian should be replaced by the Euler-Lagrange derivative of the gravitational Lagrangian as if the Riemann tensor were an independent variable, this is the form of a total derivative, Therefore, the integral can be performed without knowing the details of g(y) -we only require the conditions (21) -and the entropy reads Now, using again Eqs. (32) we can study the entropy as a function of x 0 and Q. For instance, in the limit x 0 << Q, we obtain the following approximate value, while for large x 0 we have to distinguish between the two different possibilities, where α is the parameter that we introduced in (42). The complete profile of S(x 0 ) for various values of the charge is shown in Fig. 4.
One disadvantage of this analysis is that, as we mentioned earlier, the parameter x 0 cannot be identified with the angular momentum, and therefore, the relation S(x 0 , Q) does not have a direct physical interpretation. Nevertheless, we can also study the entropy as a function of the area and of the charge, i.e, S(A, Q), and in this case the relation is meaningful since it involves three physically relevant quantities. In fact, it is interesting to check that the entropy is not only a function of the area, since it depends also on the relative amount of charge and angular momentum of the black hole. Manipulating the equations in (32), we can write the entropy (48) in the following form, where λ(A, Q) is a function obtained as a solution of the equation (52) On general grounds, for a fixed value of the area, the charge can vary from Q = 0, which would correspond to a neutral rotating black hole, to Q 2 max = A/(4π), in whose case there is no rotation and the solution is AdS 2 × S 2 . 6 It is then an interesting exercise to determine for which of these black holes of fixed area the entropy is maximal. In Fig. 5 we show the ratio S A/(4G) as a function of the charge for several fixed values of the area. First, we observe that, indeed, the entropy does not only depend on the area, but also on the charge. For Q = Q max we get S = A/(4G), since in that case the solution has no corrections. Nevertheless, when we decrease the charge leaving the area fixed -which implies that we turn on the angular momentum -the ratio between entropy and area increases. In all cases shown we see that, for a given area, a purely rotating black hole is the one that stores more information. We also observe an interesting phenomenon taking place when the area is small enough: if A < 1.91 × 4π √ µL 2 the corresponding curve becomes multivalued, indicating the existence of several black holes with same area and charge, but different entropy. This suggests the presence of a phase transition from the black hole of smaller entropy to the one of larger entropy. In that case, we see that the phase space would contain a critical point at A cr ≈ 1.91 × 4π √ µL 2 , 6 When the area is sufficiently small we obtain solutions that have Q > Qmax -see Fig. 8 -but here we focus only in the case in which Q ranges from 0 to Qmax for simplicity.  G . This picture is not completely accurate, though, because one should fix the angular momentum instead of the area in order to compare different solutions, and also because at zero temperature one cannot speak of phase transitions. Nevertheless, this result does suggest that some sort of decay could take place from one type of solution to another.

V. ADDITIONAL SOLUTIONS
In the previous section we focused on the branch of solutions that are smoothly connected to an AdS 2 × S 2 geometry, since these are particularly relevant -and the only ones that exist in Einstein gravity. However, when we solve the system of equations (32) we observe that other solutions for n(x 0 , Q) and ω(x 0 , Q) exist. A useful way of visualizing the space of solutions is to study the relation A(x 0 ) for fixed values of the charge, which we show in Fig. 6. This plot contains the curves that we showed in Fig. 3, but we see that new branches appear. In fact, for fixed values of Q and x 0 there can be up to four different solutions, which represent black holes with very different properties.

A. Branches of solutions
In the limit of x 0 → ∞, we observe that there are only two possible solutions; one which recovers the properties of extremal Kerr-Newman black holes -in particular, A → 4π(Q 2 + 2x 2 0 ) -and another one whose area tends to a constant -see Eq. (42). On the other hand, near x 0 = 0 we have in general four different solutions, which can be obtained by assuming different expansions of the parameters n and ω, as we show in the Appendix B.
One of them belongs to the AdS 2 × S 2 branch that we studied in the previous section, so we will now analyze the additional solutions.

Branch A
One possible solution of the equations (32) yields where we recall that A = 4πx 0 /ω. It is important to note that the near-horizon geometry corresponding to this choice of parameters exists for arbitrarily small values of x 0 , but not for x 0 = 0. One remarkable fact about this solution is that in the limit of x 0 → 0 the area tends to a constant value which is independent of the charge. On the other hand, the entropy can be computed using (48), and we obtain Thus, the entropy also tends to a universal constant value in the limit of vanishing x 0 , which interestingly enough corresponds to A/(2G). Observe that, for fixed values of Q and x 0 , this solution can be entropically favoured with respect to the one belonging to the AdS 2 × S2 branch. In fact, we get the following condition for small values of x 0 : However, this is not enough in order to argue that a transition from one solution to another will take place when that bound is saturated, since the angular momentum could depend differently on x 0 in both solutions, and hence, we would be comparing black holes with different conserved charges. In fact, this solution has x 0 /n ∼ 1/x 0 when x 0 → 0, which implies that the geometry departs largely from AdS 2 × S 2 , and this suggests that it actually could have a large angular momentum.
Finally, let us comment on how the black holes in this branch behave as we increase the angular momentum. Looking at Fig. 6 we observe three possibilities. If the charge is large enough, there is a maximum value of x 0 for which we can extend the branch, and at this point it merges with branch C. If the charge is smaller, the branch is connected to the solutions that have a finite area in the limit x 0 → ∞, and if it is small enough (Q < Q thr ≈ 1.13µ 1/4 L), it is connected to the Kerr-Newman branch. In other words, this implies that if we take an initial black hole with little charge but large area and angular momentum, the black hole will approach one of the solutions in this branch as it losses angular momentum, instead of an AdS 2 × S 2 geometry.

Branch B
The second additional solution has the following values of n 2 and A In this case, the area tends to zero independently of the charge when x 0 → 0. Also, unlike in the previous case, we have x 0 /n → 0, which we can interpret as a sign that the geometry is indeed slowly rotating. Now, the most interesting fact about this branch of solutions is that, even though the area vanishes in the limit of x 0 → 0, their entropy remains finite, namely Thus, the entropy per unit area in these black holes becomes arbitrarily large.

Branch C
The third and last additional solution allows for a series expansion in powers of x 1/2 0 , and the leading terms for n 2 , area and entropy read Note that, again, x 0 /n → 0, so that this solution can actually be slowly rotating, while the entropy tends to S → A/(2G).
Once the desired branch is chosen, it is possible to solve the equation (25) numerically in order to obtain the profile of g(y), as we explained previously. A comparison between these solutions is shown in Fig. 7.

B. Entropy as a function of area and charge
The preceding analysis is useful in order to characterize the space of near-horizon geometries of ECG, but it has the disadvantage that we cannot interpret x 0 as the spin parameter a. Thus, it is more meaningful to study the relation S(A, Q), which we can find exactly by using Eqs. (51) and (52). In Fig. 5 we only plotted part of this relation. The complete structure of S(A, Q) including all the solutions is quite involved and we show it as a 3-dimensional plot in Fig. 8. In obtaining this surface FIG. 8. Black hole entropy as a function of the area and charge. The thick color lines represent the different x0 → 0 limits: the red line corresponds to the AdS2 × S 2 solutions, while yellow, black and blue lines correspond to branches A, B and C, respectively. We work in units such that µL 4 = 1.
we have taken into account that the solutions of Eq. (52) must be such that n 2 > 0 and x 2 0 > 0. The red line corresponds to the AdS 2 × S 2 geometries, and interestingly these are the only ones for which S = A/(4G) -any other solution has S > A/(4G). We also represent the various x 0 → 0 limits, which correspond the yellow, black and blue curves (for branches A, B and C respectively).
As we can see, for large enough horizon area, the surface in Fig. 8 has only one branch, which recovers the Einstein gravity behaviour when A → ∞. Now, imagine that we take one of these large black holes and we start decreasing the area leaving the charge fixed -this could be interpreted as the black hole radiating away the angular momentum. 7 We find that there are two possible endpoints of this process: if the charge is large enough, then at some point we hit an AdS 2 × S 2 geometry and the black hole has radiated all the angular momentum. In order to continue evaporating it must now lose charge. On the other hand, if the charge is too small (as we saw earlier, Q < Q thr ≈ 1.13µ 1/4 L), we approach the yellow line, which corresponds to the x 0 → 0 limit of branch A, and for which A = 4π √ 3µL 2 . Interestingly, in this situation the area and the entropy of the black hole remain constant even if it loses (or gains) charge. Thus, the final product of black hole evaporation is quite different depending on which path we follow in the phase space. We also observe that for small A the surface S(A, Q) is multivalued, hence transitions or decays between solutions might occur. This illustrates that the phase space of (extremal) black hole solutions may become quite complicated in higher-derivative gravity.

VI. DISCUSSION
In this paper we have provided the first nonperturbative examples of near-horizon geometries of rotating black holes in higher-order gravity. This has been possible thanks to the special form of the equations of motion of Einsteinian cubic gravity -the density given by (2) -which can be reduced to a single second-order differential equation for one variable. Even more striking, we have been able to obtain the area and the entropy exactly in terms of the parameters of the solution, and in particular, we found the relation between black hole area, charge and entropy, S(A, Q) -see Eqs. (51) and (52). It must be noted that obtaining these quantities analytically is not possible in general higher-order theories, where the simplification of the equations that we reported does not take place. However, we do expect that there is a subset of Generalized Quasi-topological theories for which the same simplification takes place. This subset will correspond to the same type of theories admitting taub-NUT solutions that was studied in [63], where, in particular, a quartic four-dimensional density of this kind was constructed. We expect that higher-order versions of these densities exist as well, and it would be interesting to study extremal near-horizon geometries in this family of theories, thus generalizing the results presented here. In fact, we believe that the higher-order generalizations could solve some of the difficulties that we have found in our analysis and that we discuss next.

A. Large angular momentum?
Perhaps the most worrisome problem we have found is that the equation (25) seems to have no smooth solutions when the angular momentum is large compared to the charge. In particular, purely rotating (regular) black holes do not exist even in the regime where the corrections are supposed to be small. The reason, as we explained, is the vanishing of the coefficient of g in Eq. (25) at some point, which implies that the solution will not be smooth there. This issue could go away for higherorder densities, or for some appropriate combination of those, and it would be interesting to explore this possibility. On the other hand, this problem could be related to the fact that we are dealing with extremal black holes. It is known that there are certain difficulties associated with extremality (e.g., the instability of horizons [68]), and these arise explicitly in the case of higher-derivative theories -we further comment on this below. Therefore, it might happen that the problem of non-existence only affects to extremal black holes, but that (arbitrarily) near-extremal ones are fine. Despite this drawback, we believe the values found for the entropy and area of these black holes are meaningful even in the region of parameter space where no solution seems to exist. Indeed, from the point of view of EFT one should assume a perturbative expansion of the solution, and in this scheme the issue in the differential equation (25) disappears. Thus, at least the perturbative corrections to the entropy, should be meaningful in the full parameter space.

B. Multiplicity of solutions
Paradoxically enough, when Eq. (25) allows for solutions, it has many. We have seen that for fixed values of x 0 and Q, we have usually several branches of solutions with different values of the area and the entropy. But we also observed that, even when the corresponding branch has been chosen, the equation (25) can have several solutions. This is, we can have different near-horizon geometries with the same values of the charge, x 0 , area and entropy, which only differ in the shape of the horizon. Thus, in Sec. IV we only constructed numerically the solutions that are smooth deformations of AdS 2 ×S 2 geometries, but in general there are more solutions which are characterized by the same set of integration constants. In particular, the equation (36) corresponding to the limit x 0 → 0 seems to have an increasing number of solutions as Q grows. This means that there are solutions of the form AdS 2 × M 2 , where M 2 is not a sphere, but nonetheless all of these solutions have the same area and entropy. A similar situation occurs for finite x 0 . While this is an interesting phenomenon, a thorough classification of these solutions would considerable enlarge the present manuscript, and thus these additional solutions could be studied elsewhere. The degeneracy of solutions seems to be related to the sign of the higher-order coupling µ, and it would have not appeared had we taken µ < 0. The reason for taking µ > 0 is that this is required in order for asymptotically flat/AdS black holes to exist [53]. However, it is possible that for other higher-order densities the sign that allows for black holes is the same that would yield unicity of near-horizon geometries. Another relevant question is whether there exist global black hole solutions (containing an asymptotic region) of which the solutions we have constructed are the nearhorizon limit. Although it may appear shocking at first, we do not expect those solutions to exist. The reason is that the boundary problem in higher-derivative gravity is not well-posed in the presence of a degenerate horizon. This is more easily understood in the case of static, charged black holes, which allow for a simple description in ECG. Those solutions were briefly discussed in [53], where, similarly to the case here, it was shown that the equations of motion reduce to a second-order equation for one variable. Then, one has to impose a boundary condition at infinity and another one at the horizon, and this fixes the solution. But when the horizon is degenerate, the condition at the horizon turns out to fix two integration constants and it is not possible to demand the asymptotic condition. Hence, no black hole solutions exist in that case. Nevertheless, arbitrarily near-extremal ones exist, and we expect that the same behaviour will be found in the rotating case. Hence, the near-horizon geometries we have constructed make sense as a limit that non-extremal black holes can approach, but never reach. In particular, the area and entropy (and also the shape) of non-extremal black holes will tend to those found here when they approach extremality.

D. Asymptotic charges
Finally, one limitation of the near-horizon analysis is that we lose the information about the mass and the angular momentum of these black holes. We argued that the variables x 0 and n would be related, respectively, to the spin a and to the mass M but we lack a precise relation. Knowing the values of a and M would be very interesting in order to study corrections to the extremality bound and to determine the relation between the entropy and the physical charges, S(a, Q). A possible direction to achieve this goal would entail finding a generalization of Komar charge for higher-order gravities that would allow us to write the asymptotic charges as an integral over the horizon [69].

ACKNOWLEDGMENTS
We are pleased to thank Pablo Bueno, Roberto Emparan and Robie Hennigar for useful comments. PAC was mostly funded by Fundación la Caixa through a "la Caixa the opposite to the one required in order to produce inflation. However, for quartic densities (and in general, densities containing even powers of the curvature), both signs agree [61].
-Severo Ochoa" International pre-doctoral grant. In the last stages of this project PAC was supported by the KU Leuven grant "Bijzonder Onderzoeksfonds C16/16/005 -Horizons in hoge-energie fysica". DP is funded by a "Centro de Excelencia Internacional UAM/CSIC" FPI pre-doctoral grant. PAC and DP were further supported by the MINECO/FEDER, UE grant PGC2018-095205-B-I00, and by the "Centro de Excelencia Severo Ochoa" Program grant SEV-2016-0597.

Appendix A: Equations of motion
When evaluated on N (x) = 1, the gravitational tensor in (4) has the following non-vanishing components: In addition, we can relate E rr to E xx thanks to the Bianchi identity ∇ µ E µν = 0, Thus, everything is determined by the component E xx , which reads The electromagnetic energy-momentum tensor has the same structure and hence the equations of motion are reduced to E xx = T xx .

Appendix B: Solution of the Thermodynamic Quantities
Three of the four branches of solutions of the constraint equations (32) belong to the following class, where (α, β) are real parameters and we assume n 0 (Q) = 0 and ω 0 (Q) = 0. The choice (α, β) = (0, 1) corresponds to the AdS 2 × S 2 branch, while the choices (α, β) = (2, 1) and (α, β) = (1/2, 0) lead to other two solutions. The remaining branch belongs to the class The solution for all coefficients in each of the expansions can be found explicitly. In the following we exhibit the first four terms of the solutions for n 2 and ω, as well as four terms of the corresponding expansions of the area A, Wald entropy S, and relative entropy S/S BH , where S BH = A/(4G).
For the branch corresponding to AdS 2 × S 2 , determined by the coefficients (α, β) = (0, 1), For the branch determined by the coefficients (α, β) = (2, 1), For the branch determined by the coefficients (α, β) = (1/2, 0) Finally, for the class of solutions defined by (B2) π 6Q 4 − 97µL 4 6 2 3/4 3 1/4 µ 1/4 LQ 5/2 + ( Expanding g(y) as in (35) and choosing the parameter configuration of the AdS 2 × S 2 branch, we see that the solutions for all g k (y) that satisfy the boundary conditions are polynomial in y. The first terms read as follows We shall focus on the neighbourhood of solutions for which ω = 0. Such solutions are non-rotating if x 0 = 0 in such a way that lim ω→0 ω/x 0 = 0, while lim ω→0 x 0 /n = 0. On the other hand, if lim ω→0 ω/x 0 = 0 and lim ω→0 x 0 /n = 0, the solutions correspond to an ultra-spinning limit and exhibit a non-compact horizon -this only happens for Λ < 0. Let us first consider the slowly rotating case. Imposing ω = 0, one solution is x 0 = 0 and n 2 = 1 − 4ΛQ 2 − 1 − 4ΛQ 2 / −2Λ 1 − 4ΛQ 2 . Then, in a neighbourhood of this solution, ω and n 2 read, in powers of x 0 , On the other hand, performing an expansion of g(y) around x 0 = 0 as in (35), the first non-vanishing term reads Thus, we conclude that AdS 2 × S 2 is not corrected by ECG also when Λ = 0. However, it admits smooth corrections when the spin is turned on. This is analogous to the AdS 2 × S 2 branch for Λ = 0 discussed above and, in fact, taking is infinitely large, in the sense that the proper length of coordinate curves tangent to ∂ y (which extend from y = −1 to y = 1) is infinite. However, the horizon has a finite area, A = 2∆ϕ, and Wald's correction to the Bekenstein-Hawking entropy vanishes, as can be deduced from (47), because now both g(1) and g (1) are zero. Nevertheless, one can check that the profile of the solution is not going to be the same as in Einstein gravity. This solution of ECG is analogous to the super-entropic black holes of Ref. [70], in the sense that both have non-compact horizons with finite area and can be understood as an entropy-divergent limit of a rotating solution. Thus, it would be interesting to study whether this solutions do or do not respect the Isoperimetric Inequality in the context of extended black hole thermodynamics. However, further investigation in these lines is left for future work.