Spinning black holes in shift-symmetric Horndeski theory

We construct spinning black holes (BHs) in shift-symmetric Horndeski theory. This is an Einstein-scalar-Gauss-Bonnet model wherein the (real) scalar field couples linearly to the Gauss-Bonnet curvature squared combination. The BH solutions constructed are stationary, axially symmetric and asymptotically flat. They possess a non-trivial scalar field outside their regular event horizon; thus they have scalar hair. The scalar"charge"is not, however, an independent macroscopic degree of freedom. It is proportional to the Hawking temperature, as in the static limit, wherein the BHs reduce to the spherical solutions found by Sotirou and Zhou. The spinning BHs herein are found by solving non-perturbatively the field equations, numerically. We present an overview of the parameter space of the solutions together with a study of their basic geometric and phenomenological properties. These solutions are compared with the spinning BHs in the Einstein-dilaton-Gauss-Bonnet model and the Kerr BH of vacuum General Relativity. As for the former, and in contrast with the latter, there is a minimal BH size and small violations of the Kerr bound. Phenomenological differences with respect to either the former or the latter, however, are small for illustrative observables, being of the order of a few percent, at most.


Introduction
Scalar-tensor theories of gravity have attracted much attention since the pioneering example of Brans-Dicke theory [1]. The physical relevance of such models could be tested, in particular, in strong gravity systems, namely black holes (BHs). On the one hand, as it turns out, the BH solutions in Brans-Dicke theory, as well as in a large class of models where the scalar field is nonminimally coupled to the Ricci scalar, are the same as in General Relativity (GR) [2,3]. On the other hand, BHs in extended scalar-tensor models, namely those with higher curvature corrections are, generically, different from those of GR [4].
Within the class of scalar-tensor theories that possess higher curvature corrections, those including a real scalar field, φ, with a canonical kinetic term, non-minimally coupled to the Gauss-Bonnet (GB) quadratic curvature invariant, have attracted considerable interest. This is the class of Einstein-scalar-GB (EsGB) models described by the action where α is a dimensionful coupling constant and f (φ) is a dimensionless coupling function. In these models, the GB term becomes dynamical in four spacetime dimensions, and the equations of motion remain second order, which is typically not the case when higher curvature corrections are included in the action. Moreover, the GB term as a higher order correction is suggested from string theory [5].
The status of BHs in the family of models (1.2) depends on the properties of f (φ); its choice determines if φ = 0 is a consistent truncation of the equations of motion. There are two generic cases. Following the classification in [6] for a cousin model, we call models where φ = 0 is not a consistent truncation of the equations of motion class I or dilatonic-type. In this class of EsGB models φ ≡ 0 does not solve the field equations. Thus the Schwarzschild/Kerr BH is not a solution. In terms of the coupling function, this class of models obeys (from the scalar field equation (2.6) below) A representative example of coupling for this class is the standard dilatonic coupling, f (φ) = e γφ , which emerges in Kaluza-Klein theory, string theory and supergravity. In this case φ is often referred to as the dilaton field. BHs in the Einstein-dilaton-GB model were constructed in [7][8][9], where they were shown to have a qualitatively novel feature: a minimal BH size, determined by the coupling constant α. Some of these BHs are perturbatively stable [10] and aspects of their phenomenology has been considered in e.g. [11][12][13]. Models where φ = 0 is a consistent truncation are called class II or scalarised-type. In this case φ ≡ 0 solves the field equations and thus Schwarzschild and Kerr BHs are solutions of the full model. This demands that This condition holds, for instance, if one requires the model to be Z 2 -invariant under φ → −φ. The Schwarzschild/Kerr BH solution is not, in general, unique. These EsGB models may contain a second set of BH solutions, with a nontrivial scalar field profile -the scalarised BHs. Such second set of BH solutions may, or may not, continuously connect with GR BHs. Models within this class have been recenly under scrutiny in relation to BH spontaneous scalarisation -see e.g. [14][15][16][17][18].
Although f 1 is the linearisation of f 2 (the constant term is irrelevant here) these two models have qualitatively different properties. Namely, the spherical scalarised BHs with the former coupling function are unstable against perturbations; but the ones with the latter coupling function can be stable [19].
In this paper we are interested in a model of class I, the linear coupling or shift symmetric model. The coupling function is (1.5) which implies the existence of a shift symmetry: the equations of motion are invariant under the transformation with φ 0 an arbitrary constant. This invariance results from the fact that in four spacetime dimensions the GB term alone is a total divergence. BHs in the model (1.2) with (1.6) have been first discussed by Sotiriou and Zhou (SZ) [20,21]. This model falls within the Horndeski class [22], for which a no-scalar-hair theorem had been established [23]. However, the SZ solution circumvents this theorem, since one of the assumptions (finitness of a certain current) is violated. The SZ solution has a minimal size, such as the BHs in Einstein-dilaton-GB. In fact, the model (1.2) with (1.6) can be seen as a linearisation of the Einstein-dilaton-GB model, and thus one expects similar properties for the BH solutions of both models. However, as pointed out above, models with a certain coupling function and its linearisation may have different properties. It has also been argued that the SZ could emerge dynamically in a gravitational collapse scenario [24]. The goal of this paper is to construct and study the basic physical properties of the spinning generalisation of the SZ solution, which, up to now, have not been considered. Astrophysical BHs have angular momentum. Thus, considering spinning BHs is fundamental to assess the physical plausibility of any BH model. This is, however, technically more challenging than for spherical BHs, in particular in the presence of higher curvature corrections, such as the GB invariant, as described below.
This paper is organised as follows. In Section 2 we briefly discuss the equations of motion and some relevant properties of the model. In Section 3 we provide a short review of the spherical SZ solutions, as a warm up for the spinning case. In Section 4 we introduce the framework for the construction of spinning BHs, discussing the ansatz, boundary conditions, the physical quantities of interest and the numerical procedure. In Section 5 we describe the spinning BH solutions, its domain of existence, and the behaviour of different physical quantities. In Section 6 we present conclusions and remarks. Two appendices give some technical details on the construction of perturbative and extremal solutions.

The model
We consider a general EsGB model with the action (1.2). We use units such that c = 1 = 16πG. Observe that the coupling constant has physical dimension [α] ∼ [L] 2 , where L represents "length".
Varying the action (1.2) with respect to the metric tensor g µν , we obtain the Einstein field equations The effective energy-momentum tensor has two distinct components, (2. 2) The first one is due to the scalar kinetic term in (1.2) the second one is due to the scalar-GB term in (1.2), and reads where we have defined Here, ε αβρσ is the Levi-Civita tensor. The equation for the scalar field is As pointed out in the introduction, the GB term is a total divergence: where the vector P µ takes a particularly simple form [25] for a spacetime possessing a Killing vector ∂/∂t (t is the time coordinate), Thus the transformation (1.6) does not change the equations of the model. Moreover, (2.7) implies that the equation for the scalar field (2.6) can be written as As we shall see, a consequence of this relation is that the scalar 'charge' (as read off from the asymptotically leading monopolar mode) is just the Hawking temperature of BH [26].
In this work we shall be interested in stationary, axially symmetric solutions. They possess two asymptotically measured global charges: the mass M and the angular momentum J. There is also a scalar charge Q s , but it is not an independent quantity; it depends on the BH mass and angular momentum. Thus the scalar hair is of secondary type [4]. Also, note that the shift symmetry (1.6) is broken by imposing φ(∞) = 0. Horizon quantities of physical interest, on the other hand, include the Hawking temperature T H , the horizon area A H and the entropy S, whose concrete expressions are given below. Since

Spherically symmetric black holes
Before discussing the case of spinning BHs, it is of interest to review the construction and basic properties of the static, spherically symmetric BHs, the SZ solutions [20,21]. As we shall see, they contain valuable information, and share some key properties with their rotating counterparts, being easier to study since they are found by solving a set of ordinary differential equations. Moreover, a perturbative exact solution is available in the static case, which is discussed in Appendix A.1.

The equations and boundary conditions
The spherical BHs of (1.2) with (1.6) can be found using Schwarzschild-like coordinates, with a metric ansatz containing two unknown functions, where r and t are the radial and time coordinate, respectively, dΩ 2 2 is the metric on the unit round S 2 and m(r) is the Misner-Sharp mass [27], which obeys m(r) → M as r → ∞. The scalar field φ is a function of r only. The Schwarzschild BH corresponds to φ = 0, m(r) = r H /2 =constant, σ(r) = 1. One can easily verify that for α = 0 this is not a solution of the model in this work.
The advantage of this metric gauge choice is the simple form of the Einstein equations (2.1), which yield the generic relations For the considered EsGB model, the diagonal components of the effective energy-momentum tensor contain second derivatives of the metric functions N, σ. However, one can find a suitable combination of the field equations such that the functions m, σ still solve first order equations. These equations are The Einstein equations contain also a second order equation which provides a constraint, being a linear combination of (3.3) and (3.4) together with their first derivatives. The scalar field φ is a solution of a 2nd order equation in terms of N and φ only This approach leads to a good accuracy of the numerical results, and can easily be generalized for an arbitrary coupling function f (φ). The approximate form of the solutions valid for large-r reads (3.6) in terms of mass M and a scalar "charge" Q s . Close to the event horizon, located at r = r H , the solutions possess an approximate expression as a power series in r − r H , with where while φ 2 is a complicated function of φ 1 , r H and α. The Hawking temperature, horizon area and entropy of the solutions, as computed from the formalism in the next Section, are given by The field equations imply that the first derivative of the scalar field, φ 1 , is a solution of the quadratic equation which implies the following condition for the existence of a real root This requirement translates into the following coordinate independent condition between the horizon size and the coupling constant α We remark that A H = 4πr 2 H for the metric ansatz employed here. Thus, for a theory with a given value of the input parameter α > 0, the BHs are not smoothly connected with the Minkowski vacuum. There is minimal horizon size and a mass gap [20,21], just as for BHs in the Einsteindilaton-GB model [7][8][9].

The solutions
The parameter space of solutions can be scanned by starting with the Schwarzschild BH (α = 0) and increasing the value of α for fixed r H . When appropriately scaled, they form a line, starting from the smooth GR limit and ending at a critical solution where the condition (3.12) is violated, and where the maximal value of the ratio α/M 2 (around 0.32534) is achived. Once the critical configuration is reached, the solutions cease to exist in the parameter space. Physically this means that the EsGB BHs have a minimal size and mass, for given α. A possible interpretation is that the GB term provides a repulsive contribution, becoming overwhelming for sufficiently small BHs, thus preventing the existence of an event horizon. The full set of static solutions will be shown below in Fig. 3 (the blue dotted line with j = 0) as a function of the dimensionless parameter α/M 2 .
As discussed in Appendix A.1, a simple perturbative solution can be found as a power series in the parameter The results in Appendix A.1 imply the following expressions Interestingly, all corrections to the reduced temperature t H are positive. That is, for the same mass, the shift symmetric Hordenski BH is 'hotter'. For the other quantities, no clear generic pattern emerges.
We have found that the perturbative solution provides a very good approximation to the numerical results. This follows from the smallness of the parameter β. In fact, condition (3.11) implies β max 0.102062. As such, the contribution of the higher order terms in β quickly becomes irrelevant.

Ansatz and boundary conditions
To obtain stationary and axi-symmetric BH spacetimes, possessing two commuting Killing vector fields, ξ and η, we use a coordinate system adapted to these symmetries. Then ξ = ∂ t , η = ∂ ϕ , and we consider a metric ansatz which has been employed in the past for the study of Kerr BHs with scalar hair [28]. In terms of the spheroidal coordinates r, θ and ϕ (with t the time coordinate), the metric line element reads: where the metric functions F i , W , as well as the scalar field φ, depend on r, θ only and r H > 0 is an input parameter again describing the location of the event horizon. The coordinates θ, ϕ and t possess the usual range, while r H r < ∞. The vacuum Kerr BH can be written in this form, the corresponding expressions of F 0 , F 1 , F 2 and W being displayed in Appendix A of [29]. Finding BH solutions with this ansatz requires defining boundary behaviours. We have made the following choices. For the solutions to approach at spatial infinity (r → ∞) a Minkowski spacetime we require lim Since the scalar field is massless, one can construct an approximate solution of the field equations compatible with these asymptotics as a power series in 1/r. The leading order terms of such an expansion are: where c t , c ϕ and Q s are constant parameters to be fixed by the numerics. Axial symmetry, together with regularity at the axis impose the following boundary conditions on the symmetry axis, i.e. at θ = 0, π: As before, an approximate expansion of the solution compatible with these boundary conditions can be constructed; as an illustration, at θ = 0 one finds where F a = {F 0 , F 1 , F 2 , W ; φ}. The essential data, which is fixed by the numerics, is encoded in the functions F a0 = {F i0 , W 0 , φ 0 }. Moreover, the absence of conical singularities implies also that F 1 = F 2 on the symmetry axis. Focusing on BHs with parity reflection symmetry, we need to consider the solutions only for 0 θ π/2. Then, the functions F i , W and φ satisfy the following boundary conditions on the equatorial plane (θ = π/2) For the metric ansatz (4.1), the event horizon is located at a surface with constant radial variable, r = r H > 0. By introducing a new radial coordinate the horizon boundary conditions and numerical treatment of the problem simplify. These boundary conditions are where Ω H is the horizon angular velocity, and the Killing vector χ = ξ + Ω H η is orthogonal and null on the horizon. These conditions are consistent with the near horizon solution where the essential functions are F i0 (also F 0 r H = F 1 r H ).

Quantities of interest and a Smarr relation
Many quantities of interest are encoded in the metric functions at the horizon or at infinity. Considering first horizon quantities. The Hawking temperature is The horizon angular velocity Ω H is fixed by the horizon value of the metric function W , The total (ADM) mass M and angular momentum J of the BHs are read off from the asymptotics of g tt and g ϕt , These global quantities can be split into the horizon and bulk contributions -see, e.g., [30]. These are, respectively M H and J H , computed as a Komar integrals on the horizon, and M φ and J φ , computed as volume integrals of the appropriate effective energy-momentum tensor components: where Σ is a spacelike surface, bounded by the 2-sphere at infinity S 2 ∞ and the spatial section of the horizon H. M φ and J φ encode the contribution of the effective "matter" distribution to the total mass and angular momentum. For Kerr BHs, M = M H and J = J H ; this is not so for EsGB BHs.
= 0, only the GB part of the effective energy-momentum tensor (2.2) contributes to the energy and angular momentum "matter" densities.
The solutions can be shown to obey the Smarr-type law where S is the entropy as computed from Wald's formula [31], and R is the Ricci scalar of the induced horizon metric h. In the Smarr-type law, M s is a contribution of the scalar field which can also be expressed as an integral of φR 2 GB term. Also, by integrating (2.9) over an hypersurface bounded by the event horizon and the sphere at infinity one can prove the following interesting relation This proportionality between the scalar charge and the Hawking temperature is a unique feature of the shift symmetric EsGB model, see the discussion in [26]. The EsGB BHs satisfy also the first law

The numerical approach
In our approach, the field equations reduce to a set of five coupled non-linear elliptic partial differential equations for the functions F a = (F 0 , F 1 , F 2 , W ; φ), which are found by plugging the ansatz (4.1) together with φ = φ(r, θ) into the field eqs. (2.1), (2.6). They consist of the Klein-Gordon equation (2.6) together with suitable combinations of the Eintein equations (2.1) The explicit form of the equations solved in practice is too complicated to display here; each equation containing around 250 independent terms. Also, the remaining equations E r θ = 0 and E r r − E θ θ = 0 are not solved directly, they yielding two constraints which are monitored in numerics. Typically they are satisfied at the level of the overall numerical accuracy. We remark that one can verify that the remaining equations vanish identically, E ϕ r = E t r = E ϕ θ = E t θ = 0, the circularity condition being satisfied. As such, the employed ansatz is consistent, a fact which is not a priori guaranteed (see [32] for a discussion in an Einstein-scalar field model which leads to a non-circular metric form).
Our numerical treatment can be summarised as follows. We restrict the domain of integration to the region outside the horizon. Then, the first step is to introduce the new radial variablē x = x/(1 + x) which maps the semi-infinite region [0, ∞) to the finite region [0, 1], where x is given by (4.7) and r is the radial variable in the line element (4.1). Next, the equations for F a are discretised on a grid inx and θ. Most of the results in this work have been found for an equidistant grid with 300 × 40 points. The grid covers the integration region 0 x 1 and 0 θ π/2.
The equations for F a have been solved subject to the boundary conditions introduced above. All numerical calculations are performed by using a professional package [33], which employs a Newton-Raphson method. This code uses the finite difference method, providing also an error estimate for each unknown function. For the solutions in this work, the maximal numerical error for the functions is estimated to be on the order of 10 −3 . The Smarr relation (4.15) provides a further test of the numerical accuracy, leading to error estimates of the same order.
In our numerical scheme, there are three input parameters: i) the event horizon radius r H ; ii) the event horizon angular velocity Ω H in the metric ansatz (4.1) and iii) the coupling constant α in the action (1.2). The quantities of interest are computed from the numerical output. For example, the mass M , and the angular momentum J are extracted from the asymptotic expressions (4.12), while the Hawking temperature, the entropy and the horizon area are obtained from the event horizon data.
The results reported in this work are obtained from around twenty thousand solution points. For all these BHs we have monitored the Ricci and the Kretschmann scalars, and, at the level of the numerical accuracy, we have not observed any sign of a singular behaviour on and outside the horizon (see, however, the discussion below on the limiting solutions).

General properties and limiting behaviour
In an approach based on the Newton-Raphson method a good initial guess for the profile of the various functions is an essential condition for a successful implementation. The spinning solutions in this work can be constructed by using two different routes. In the first approach, one uses the profile of a Kerr BH with given r H , Ω H as an initial guess for EsGB solutions 1 with a small value of the ratio α/r 2 H . The iterations converge and, repeating the procedure, one obtains in this way solutions with large α. In the second approach, one starts instead with spherically symmetric solutions of EsGB, either obtained numerically or from the perturbative expansion. These can also be studied within the ansatz (4.1), with W = 0, F i being functions of r only and with F 1 = F 2 . Then, starting with an EsGB spherical BH with a given r H and α = 0, rotation is introduced by introducing and slowly increasing Ω H .
For all solutions we have found, the metric functions F a , together with their first and second derivatives with respect to both r and θ have smooth profiles. This leads to finite curvature invariants on the full domain of integration, in particular at the event horizon. The shape of the metric functions F 0 , F 1 , F 2 and W is similar to those in the α = 0 case. The maximal deviation from the Einstein gravity profiles (with the same input parameters r H , Ω H ) is near the horizon. At the same time, the scalar field may possess a complicated angular dependence, in particular for fast spinning configurations.
The profile functions of a typical solution are exhibited in Figure 1. The insets show the same curves for Kerr with the same r H , Ω H , for comparison. The Ricci and the Kretschmann scalars, R and K ≡ R αβµν R αβµν , together with the components T t t and T t ϕ of the effective energy-momentum tensor are shown in Figure 2. In these plots, the corresponding functions are shown in terms of the (inverse) radial variable r for three different values of the angular coordinate θ. One observes, for instance, that g tt becomes positive along the equator, near the horizon, thus manifesting the existence of an ergo-region (see next subsection). One also notices that both R and K stay finite everywhere, in particular at the horizon. From the components of the effective energy-momentum tensor one observes, in particular, that −T t t < 0 for a region in the vicinity of the symmetry axis, manifesting a breakdown of the weak energy condition for the effective energy-momentum tensor.
Returning to the construction of the solutions, we have noticed the existence of a critical set of input parameters for which the numerical process fails to converge. Neither a singular behaviour nor a deterioration of the numerical accuracy in the vicinity of this set was observed. An explanation for this behaviour, similar to that justifying the critical configurations found in the static case, is based on the analysis of the field equations in the vicinity of the event horizon. After some algebra, one finds that the second order term φ 2 (θ) in the expansion of the scalar field φ(x, θ) = φ 0 (θ) + φ 2 (θ)x 2 + . . . is a solution of a quadratic equation, where the coefficients a, b, c depend on the values of F i , W and their derivatives at the horizon. Then, a real solution to the above equation exists only if ∆ = b 2 − 4ac > 0. In practice, we have monitored this discriminant and observed that the numerical process fails to converge 2 when ∆ takes small values close to zero at θ = 0, π. As in the spherically symmetric case, we have found no evidence for the emergence of a secondary branch of solutions in the vicinity of the critical solutions.
A different limiting behaviour is found when varying the value of the horizon velocity Ω H for fixed (r H , α). As for the vacuum Kerr family, following this method one finds two branches of solutions, which join for a maximal value of Ω H . The first branch emerges from the corresponding static configuration. The second branch, on the other hand, ends, as for α = 0, at extremal configurations. These have vanishing Hawking temperature and nonvanishing global charges, horizon area and entropy. We must emphasise, however, that only near extremal solutions, as opposed ot exactly extremal BHs, can be constructed within the framework proposed in this work. As such, the results for the extremal solutions reported here result from extrapolating the data found in the near-extremal case. Moreover, unlike the extremal vacuum Kerr BH which yields a perfectly regular geometry [49], the extremal EsGB solutions appear to not be regular, with the Ricci scalar tending to diverge at the poles of the horizon. A partial understanding of this behaviour is given in Appendix B, based on a perturbative construction of the near-horizon configurations.

The domain of existence
Let us now address the domain of existence of the EsGB solutions. There are two fundamental scales, the coupling constant α, and the BH mass of the solutions M . In what follows we display various quantities of interest as a function of the dimensionless coupling constant α/M 2 . This parameter measures the impact of non-GR features, due to the GB contribution. The analysis is also performed in terms of the dimensionless angular momentum j = J/M 2 . This parameter measures the impact of non-staticity. The link between these two quantities is provided by the Figure 3, where we plot the domain of existence (shaded blue region) in a j vs. α/M 2 plot. Therein, all data points which were found numerically are also explicitly shown. The blue shaded region is the extrapolation of these points into the continuum. The figure shows that the domain of existence is delimited by: • the set of static BHs (j = 0, blue dotted line); • the set of extremal BHs (black dotted line); • the set of critical solutions (green line); • the set of GR solutions -the Kerr/Schwarzschild BHs (α/M 2 = 0, red line). Two comments on Figure 3. First, the Kerr bound j 1 is violated for spinning EsGB BHs in a small region of the domain of existence close the extremal set. However, this violation is rather small, with j (max) ∼ 1.013 for all (accurate enough) solutions studied so far. Second, along j fixed lines, the critical solution is attained at a smaller α/M 2 as j is increased. A possible interpretation is that both the GB contribution and the spin are repulsive effects. Thus, in the presence of rotation, BHs cease to exist for a smaller GB contribution.
In Figure 4  Let us comment on some features resulting from Figure 4. For fixed j, the BH area decreases as α/M 2 increases; but the corresponding reduced BH entropy increases. This provides a clear example how BH entropy deviates from the Hawking-Bekenstein formula in this modified gravity: when the GB contribution becomes larger, the BH becomes smaller but it carries more entropy (for fixed j). On the other hand, fixing the EsGB dimensionless coupling constant α/M 2 , both the reduced area and the reduced entropy decrease as j increases. Thus, for any fixed EsGB model, spin reduces the size and the entropy of BHs. The BH temperature, on the other hand, increases with α/M 2 for fixed j and decreases with j for fixed α/M 2 .  For the Kerr BH, this surface has a spherical topology and touches the horizon at the poles. As discussed in [34], the ergoregion can be more complicated for other models, notably for BHs with synchronised scalar hair, with the possible existence of an additional S 1 × S 1 ergo-surface (ergotorus) -see also [35]. We have found that this is not the case for EsGB BHs, where all solutions are Kerr-like in the sense they possess a single topologically S 2 ergosurface. Let us now consider the horizon geometry. Similarly to the GR Kerr solution, EsGB BHs have an event horizon of spherical topology. The metric of a spatial cross-section of the horizon is

Ergoregion and horizon properties
Geometrically, however, the horizon is a squashed, rather than round, sphere. This is shown by computing the horizon circumference along the equator, L e , and along the poles, L p : The ratio of these two circumferences define the sphericity [36] s ≡ L e L p . (5.5) In Figure 5 (left panel) the sphericity is shown as a function of the dimensionless coupling constant α/M 2 . An interesting feature there is that s can exceed the maximal GR value for a set of EsGB solutions close to extremality. Roughly, the EsGB can become more oblate than Kerr. Also, as expected, the squashing of the horizon produced by the rotation is such that s is always larger than unity. That is, the solutions are always deformes towards oblatness, rather than prolatness. Another physical quantity of interest is the horizon linear velocity v H [36][37][38]. v H measures how fast the null geodesics generators of the horizon rotate relatively to a static observer at spatial infinity. It is defined as the product between the perimetral radius of the circumference located at the equator, R e ≡ L e /2π, and the horizon angular velocity Ω H , As seen in Figure 5 (right panel), all studied EGBs solutions have v H < 1, just like for Kerr, and despite the (small) violations of the Kerr bound. Thus, the null geodesics generators of the horizon rotate relatively to the asymptotic observer at subluminal speeds. Further insight into the horizon geometry is obtained by considering the isometric embedding of the spatial sections of the horizon in an Euclidean 3-space E 3 . A well-known feature of the Kerr horizon geometry is that for a dimensionless spin j > √ 3/2 ≡ j (S) (dubbed Smarr point) the Gaussian curvature of the horizon becomes negative in a vicinity of the poles [39]. In this regime, an isometric embedding of the Kerr horizon geometry in E 3 is no longer possible. As expected, this feature also occurs also for the solutions in this work, even though the position of the Smarr point now depends on the value of the dimensionaless coupling constant α/M 2 . Following [36,37], the collection of Smarr points as α/M 2 is varied is dubbed the Smarr line. Figure 5 displays also the position of the Smarr line as a function of α/M 2 . One observes that, as for the Kerr limit, an isometric embedding of the horizon geometry in E 3 is possible only up to a maximal value of s and v H . Also, notice that both the sphericity s and v H are not constant along the Smarr line and slighly larger values of both these quantities are allowed for embeddable BHs when α/M 2 is increased.

Orbital Frequency at the ISCO and Light Rings
A phenomenologically relevant aspect of any BH concerns the angular frequency at both the innermost stable circular orbit (ISCO) and the light ring (LR). The former is associated to a cut-off frequency of the emitted synchrotron radiation generated from accelerated charges in accretion disks. The latter is related to the real part of the frequency of BH quasi-normal modes [40]. The LRs are also key in determining the BH shadow [41]. Following a standard method, one finds that the angular frequency of a test particle with energy, E, and angular momentum, L, on the equatorial plane, θ = π/2, is, The radial coordinate, r, of such particle obeys the equation, where the 'dot' denotes derivative with respect to an affine parameter. is a constant with = 0 for massless test particles and = −1 for the massive test particles. The former are relevant for the LRs and the latter for the ISCO.
In the case of massive test particles, circular orbits require that both the potential V (r) and its derivative vanish, V (r) = V (r) = 0. This yields two algebraic equations for E and L, which can be solved analytically. These have two distinct pairs of solutions, (E + , L + ) and (E − , L − ), corresponding, respectively, to co-rotating and counter-rotating orbits. It is then possible to assess the stability of the circular orbits by computing the second derivative of the potential. The ISCO will correspond to the orbit in which the test particle has energy and angular momentum that solves V (r) = V (r) = 0 and the radial coordinate that solves V (r) = 0. Having obtained the energy, angular momentum and radial coordinate of the ISCO, the corresponding angular frequency is computed using (5.7).
In Fig. 6, we present the ratio between the angular frequency at the ISCO between EsGB BHs and Kerr BHs, for both co-rotating, ∆ω co ISCO and counter-rotating orbits, ∆ω counter ISCO , fixing j, as a function of the reduced coupling constant, α/M 2 : . (5.9) Several illustrative values of j are exhibited. For both the co-rotating and counter-rotating cases, by definition, the ratio converges to unity in the Kerr limit. For all fixed j and for both co and counter-rotating orbits, the ratio diverges away, monotonically, from unity as α/M 2 increases. How the ratio goes away from unity depends, however, on j and on the direction of the orbital motion.
For j = 0 the distinction between co and counter rotating orbits is meaningless. The ratio grows away from unity as α/M 2 increases -solid blue line in Fig. 6. Naively, this is related to the fact that the static BH size decreases with increasing α/M 2 , making the ISCO also decrease and hence its frequency increase. Introducing j raises the degeneracy between co and counter rotating orbits. For co-rotating (counter-rotating) orbits and small j = 0, the ratio is always larger (smaller) than that for the static BHs (j = 0) -dotted lines in Fig. 6 (left and right panels). One may interpret these behaviours as a consequence of frame dragging, which enhances (damps) motion along co-rotating (counter-rotating) orbits. In the counter-rotating case this trend remains for large j -dashed lines in Fig. 6 (right panel). In the co-rotating case, however, an unexpected behaviour emerges. For sufficiently large j the ratio stops being enhanced with respect to the static case, and eventually becomes suppressed with respect to it -dashed lines in Fig. 6 (left panel).
A possible explanation for this unexpected behaviour is found by studying the angular velocity of the horizon, Ω H . This quantity is a better measure of dragging effects than the spacetime angular momentum. Indeed, the fact that a BH has a large j does not imply that it has a large horizon angular velocity. 3 Let us then consider the reduced horizon angular velocity, ω H ≡ Ω H M , and its difference beween EsGB and Kerr BHs with the same j, defined as: This quantity is plotted against the reduced angular momentum j in Fig. 7. One observes that, for small enough fixed j, the EsGB BHs have larger ω H than Kerr ones. This support the thesis that dragging effects are stronger and should enhance the angular frequency at the ISCO. However, after a given spin j, the EsGB BHs have smaller ω H than Kerr BHs. That is, albeit having a larger spacetime angular momentum, large j EsGB BHs spin more slowly, and thus source weaker frame dragging, than Kerr BHs. Qualitatively, at least, this provides an explanation for the behaviour observed in Fig. 6 (left panel). Quantitatively, for co-rotating orbits, the maximal deviation from Kerr is ∆ω co ISCO ∼ 8% and occurs for j ∼ 0.5 and the maximal value of α/M 2 . For counter-rotating orbits, on the other hand, the ratio is maximised, for any α/M 2 , by the static case.
In the case of massless particles, a similar analysis can be done. Now, solving V (r) = 0, we obtain an algebraic equation for the impact parameter, b p = L/E, which yields two distinct solutions b + p and b − p corresponding to co-rotating and counter-rotating orbits, respectively. Using this result, and solving V (r) = 0, yields the radial coordinate of the LR. Having computed the impact parameter and the radial coordinate of the LRs, one can again compute their angular frequency, using (5.7). Fig. 8 shows the ratio between the angular frequency at the LR of EsGB BHs and Kerr BHs, for both co-rotating and counter-rotating orbits, ∆ω counter LR , defined in an analogous way to (5.9), with different values of spin, j, as a function of the reduced coupling constant, α/M 2 . The overall behaviour is very similar to the one discussed above for the ISCO frequency. The main difference for the LR case is that the maximal deviation for both types of orbits is smaller than the corresponding orbits at the ISCO.
to introduce this notion in BH physics.

Conclusions and remarks
In this work we have constructed the spinning generalisations of the static BHs in the shift symmetric Hordenski model. This is a family of asymptotically flat, stationary, axially symmetric BHs, that are non-singular on and outside an event horizon. The domain of existence of these solutions is naturally described by two dimensionless parameters: the dimensionless coupling constant of the model, α/M 2 , and the dimensioness spin of the BHs, j = J/M 2 . Then, the domain of existence is bounded by four special limiting behaviours: the GR limit (when α = 0), the static limit (when j = 0), the extremal limit, when the surface gravity of the solutions vanishes, and a critical set of solutions for which a horizon ceases to exist. This last boundary has an important implication. For non-zero α it means there is minimum mass (and hence) size for BHs. Thus there is a mass gap with respect to the Minkowski vacuum, which is also a solution of the theory. This non-GR property also occurs for the Einstein-dilaton-GB model discussed e.g. in [8,9]. Other properties of the BHs we have constructed and analysed in this paper also parallel the solutions found in the Einstein-dilaton-GB model. This similarity of properties was antecipated by the observation made in the introduction: the linearisation of the action of Einstein-dilaton-GB model reduces to (1.2) in the limit of small φ, i.e. e φ 1 + φ, by virtue of (2.7). Since the scalar field takes rather small values for typical Einstein-dilaton-GB BHs, the shift symmetric EsGB BHs with the same input parameters provide a reasonable approximation -see, for instance, the bottom left panel of Figure 1 for the scalar field magnitude of a typical solution. Thus, the domain of existence of the Einstein-dilaton-GB and EsGB BHs are indeed quite similar, as confirmed by the results in this work.
Yet, there are both qualitative and quantitative differences between the two models. An intriguing property of the model we have focused on, that does not occur for the Einstein-dilaton-GB model, is the scalar charge-temperature relation (4.18). In fact, also the Smarr law is different in both models. Quantitatively, the correspondence between the two models holds only for small enough values of α/M 2 and j. For example, the critical value of the ratio α/M 2 is 0.3253 for the spherically symmetric solutions in this work (being fixed by an algebraic condition between the horizon size and the coupling constant α, Eq. (3.12)) and 0.1728 for Einstein-dilaton-GB BHs (in which case the generalization of (3.12) includes, as well, a dependence on the value of the scalar field at the horizon, see e.g. Ref. [7]). Moreover, a specific feature of the Einstein-dilaton-GB model is the occurrance, near the critical configuration of a small secondary branch of BH solutions [42][43][44]. Along this branch, the mass increases with decreasing horizon radius. This secondary branch appears to be absent in the EsGB case.
Finally, let us remark that the way the SZ solutions circumvent the no-scalar-hair theorem also applies to the model herein [23]. This occurs by violating the assumption that the current associated to the shift-symmetry should be finite at the horizon. For the static SZ BHs this current diverges on the horizon. This, however, does not induce any physical pathologies. We have checked that this current (squared) diverges at the horizon also in the spinning BHs reported in this work.
The choice (A.1) leads to a particularly simple structure of the equations for the functions {h k (r), σ k (r), φ k (r)}, which can easily be solved to an arbitrary order. We have done it up to k = 12. These functions are polynomials in x = r H /r, the expression of the first few terms being h 0 (r) = 1, h 1 (r) = 0, h 2 (r) = − 49 5 x − 29 5 Unfortunately, no general pattern can be found and the coefficients of the terms in the polynomial expressions of {h k (r), σ k (r), φ k (r)} become increasingly complicated, with higher powers of x = r H /r. The corresponding expression of the mass function m(r) follows directly from (A.1) (we recall that N = 1 − 2m(r)/r). While m(r) is strictly positive, its derivative becomes negative in a region close to the event horizon, the lowest order term being

B. The attractors and the issue of extremal solutions
The numerical results suggest that, unlike the extremal Kerr solution, the extremal EsGB solutions are not regular. Evidence for this conjecture is obtained as follows. Instead of solving the full bulk equations searching for extremal solutions, one tackles the construction of the corresponding nearhorizon configurations. In this case, one has to solve a co-dimension one problem (the radial dependence being factorized), whose solutions are easier to study. Since this problem was already considered in a more general context [48] (see also the corresponding Einstein-dilaton-GB computation in [8]), in what follows we shall review the basic results only. The idea to consider a construction of the near-horizon limit of the extremal rotating BH as a power series in α. The background solution is taken to be the vacuum near horizon extremal Kerr (NHEK) solution in pure Einstein gravity [49]. As we shall see, the α 2 -corrections to this solution are singular and destroy its smoothness.
In the next step, we find the expression of φ 1 (u) by solving the eq. (2.6) in the NHEK background (B.7) such that φ 1 (u) is divergent at u = 1 only.
Then, with the above expressions, one can prove the existence of a singularity at the poles of the horizon, with the Ricci scalar diverging at θ = 0 (i.e. u = 1) Finally, let us mention that the above perturbative result does not exclude the existence of regular solutions for α large enough. Thus we have also attempted to solve the field equations of the model within a nonperturbative approach, by solving a boundary value problem. The imposed boundary conditions assure the regularity of the configurations at u = ±1. However, no such solutions could be found.