Exploring the Small Mass Limit of Stationary Black Holes in Theories with Gauss-Bonnet Terms

In this work we examine the small mass limit of black holes, with and without spin, in theories where a scalar field is non-minimally coupled to a Gauss-Bonnet term. First, we provide an analytical example for a theory where a static closed-form solution with a small mass limit is known, and later use analytical and numerical techniques to explore this limit in standard scalar-Gauss-Bonnet theories with dilatonic, linear and quadratic-exponential couplings. In most cases studied here, we find an inner singularity that overlaps with the event horizon of the static black hole as the small mass limit is reached. Moreover, since solutions in this limit possess a non-vanishing Hawking temperature, a naked singularity is expected to be reached through evaporation, raising questions concerning the consistency of these theories altogether. On the other hand, we provide for the first time in this context an example of a coupling where the small mass limit is never reached, thus preferred from the point of view of cosmic censorship. Finally, we consider black holes with spin and numerically investigate how this changes the picture, using these to place the tightest upper bounds to date on the coupling constant for the dilatonic and linear theories, with $\sqrt{\overline{\alpha}}<1$ km.

In this work we examine the small mass limit of black holes, with and without spin, in theories where a scalar field is non-minimally coupled to a Gauss-Bonnet term. First, we provide an analytical example for a theory where a static closed-form solution with a small mass limit is known, and later use analytical and numerical techniques to explore this limit in standard scalar-Gauss-Bonnet theories with dilatonic, linear and quadratic-exponential couplings. In most cases studied here, we find an inner singularity that overlaps with the event horizon of the static black hole as the small mass limit is reached. Moreover, since solutions in this limit possess a non-vanishing Hawking temperature, a naked singularity is expected to be reached through evaporation, raising questions concerning the consistency of these theories altogether. On the other hand, we provide for the first time in this context an example of a coupling where the small mass limit is never reached, thus preferred from the point of view of cosmic censorship. Finally, we consider black holes with spin and numerically investigate how this changes the picture, using these to place the tightest upper bounds to date on the coupling constant for the dilatonic and linear theories, with √ α < 1 km.

I. INTRODUCTION
Prominent examples of alternative theories of gravity to General Relativity (GR) are Einstein-scalar-Gauss-Bonnet (EsGB) theories, where a new fundamental scalar is non-minimally coupled to the Gauss-Bonnet (GB) term Such models belong to the Horndeski class of theories [1,2], and in the simplest case their action takes the form where φ is a real scalar field, ξ (φ) is the (non-minimal) coupling function, and α the GB coupling constant with dimensions of length squared. EsGB theories are of wide interest and have been the subject of many works in recent years (see e.g. ). From a fundamental physics perspective they arise, for example, as the low-energy limit of some string theories [29][30][31][32] where the scalar, the dilaton, couples exponentially to the GB term, ξ (φ) ∼ e γφ [18-20, 28, 33]. From a more phenomenological point of view, they are one of a variety of theories [3-16, 18-20, 33-36] that evade the no hair conjecture [37][38][39][40][41][42][43][44] (see [45] for a review), raising the exciting possibility they can be constrained by black hole (BH) physics in the strong-curvature regime. Moreover, in some EsGB theories a dynamical mechanism called spontaneous scalarization [6-8, 12, 14, 15, 33-36, 46, 47] can occur, such that deviations from GR occur only in the strong-curvature regime. Constraints can be theoretical, such as self-consistency, or observational in nature. Indeed, only recently has GR begun to be tested observationally in the strong-field regime [48], with the dawn of gravitational wave (GW) astronomy [49][50][51].
In a striking difference to GR, some EsGB black holes (e.g. with linear and dilatonic couplings) are known to possess a minimum mass solution whose Hawking temperature is finite and non-vanishing, naturally raising the question of what is the fate of black holes in EsGB theories (see e.g. [4,18,28,[52][53][54][55][56]). This conundrum is often overlooked, but was recently explored in Refs. [57,58], where non-linear numerical simulations of evaporating EsGB dilatonic black holes were performed, supporting the idea that the end-point of Hawking evaporation is likely a naked singularity, violating weak cosmic censorship. This is a rather concerning scenario, raising questions about the consistency of EsGB models.
In this work our main purpose is to explore further the small mass limit of black holes for several couplings in EsGB theories, including those allowing for spontaneous scalarization, where a detailed analysis of the small mass limit is so far lacking. Our aim is to investigate self-consistency and observational constraints imposed on Gauss-Bonnet theories and their coupling dependence. After providing a novel example of a closed-form solution with a small mass limit, we will complement previous studies [4,18,28,[52][53][54][56][57][58] on EsGB black hole solutions with a thorough analytical and numerical exploration of the domain of existence of solutions and their inner structure, linking the existence of an inner singularity to the repulsive effects originating from the Gauss-Bonnet term, and to the structure of the field equations. Using analytical arguments, this singularity will be shown to overlap with the event horizon in the small mass limit. We will also provide for the first time, in this context, an example of a coupling function (quadratic-exponential coupling in Eq. (22) below, with β above a certain value) where the small mass limit is never reached, showing that a minimum size solution is not a generic feature of EsGB theories 1 . Finally, we construct stationary black hole solutions by numerically solving the field equations in axi-symmetry, with the aim of exploring the small mass limit once spin is considered, and finish by imposing the tightest constraints on the coupling constant, to date, on the dilatonic and linear theories.
The rest of the paper is organized as follows. First, in Section II, we introduce a theory which, unlike Eq. (2), admits a known analytical example of a static black hole with a small mass limit. This allows us to explore key features with the advantage of an exact solution. In Section III, we discuss the form of coupling functions in standard Gauss-Bonnet theories and their corresponding different phenomenologies. Then in Section IV, we explore the small mass limit of static black hole solutions for these theories, and later impose upper bounds on the coupling constant α in Section V. Finally, in Section VI, we consider how spin changes the picture. We conclude in Section VII.

II. EXPLORING THE SMALL MASS LIMIT: AN ANALYTICAL EXAMPLE
No analytic closed-form black hole solution is known to EsGB models described by the action (2). In the existing literature, therefore, the study of critical solutions in such theories has been performed by resorting to numerical techniques. In this section, we will explore an illustrative example of a related theory with known closed-form black hole solutions.
The theory is known as gravity with a generalized conformal scalar field, and was first derived in Ref. [59] by imposing conformal invariance on the equation of motion of the scalar field. Its action is given by which contains the terms present in the EsGB model of Eq. (2) with a logarithmic coupling, as well as other non-trivial interactions including a conformal coupling to gravity and a conformally invariant quartic self-coupling. The above theory is intimately connected with the scalar-tensor formulations of so called 4D-Einstein-Gauss-Bonnet gravity [60][61][62][63][64] (see Ref. [65] for a review), and belongs to the Horndeski class with functions specified by [1,2] Note that for the purpose of presentation, we have flipped the sign of α relative to the presentation of Ref. [59], and that of typical 4D-Einstein-Gauss-Bonnet studies. Among the many interesting features of Eq. (3), one that stands out is that a special combination of the field equations decouples from the scalar field, imposing a proportionality condition between the Ricci and GB scalars This allows for an easy search of closed form solutions. One known closed-form black hole solution to the above theory (with λ −1 = 6α) is given by [59] 2 with where M is the ADM mass of the black hole, and we take α to be positive. The scalar field can be seen to be regular on and outside the event horizon located at Analysing the Ricci scalar of the solution (7) we observe that revealing the existence of two physical singularities, one located at r = 0, and a finite radius singularity located at the point where the quantity inside the square-root in Eq. (7) for the function f (r) vanishes To ensure physical behaviour of the solution we require that i) the singularity located at r = r s is hidden behind the event horizon (r s < r H ); ii) the metric functions and the scalar field in Eq. (7) are real. Under these requirements, it can be shown using Eqs. (7), (8), and (10) that the following condition must hold or, in terms of r H , In other words there is a minimum mass M min = √ α 2 √ 2 (or equivalently, a minimum horizon radius r min H = √ 2α), below which solutions can no longer be described by black holes. For an object with r H = r min H , r s and r H overlap, as can be observed in Fig. 1. A possible physical interpretation for the minimum mass solution is related to the repulsive effect of the Gauss-Bonnet term on the solutions. Examining the components of the effective stress-energy tensor we get where ρ and p r are interpreted as the effective energy density and radial pressure, respectively. 3 The effective radial pressure is positive everywhere (repulsive), and diverges at r = r s . Indeed, for the minimum mass solution, the repulsive effects of the Gauss-Bonnet term dominate over the standard attractive ones, impeding the existence of a regular horizon. An interesting remark can be made regarding the Hawking temperature of the black hole solution given by This is that not only is the temperature non-zero as r H → r min H , but in fact diverges! Therefore, as the small mass limit is approached, evaporation will not halt and the black hole will continue to lose its mass at a rate [67] where G m (ω) are the graybody factors for modes with frequency ω, angular dependence ( , m), and the plus/minus sign is related to the emission of fermions/bosons. This intriguing feature casts doubts on the endpoint of Hawking evaporation, with the above calculations suggesting that a naked singularity is a strong endpoint candidate. To summarize, there are a few lessons to be learned from the example. First, it appears that theories with highercurvature terms are susceptible to the existence of a finite radius singularity located at r s > 0. The existence of this singularity is, mathematically, intimately tied to terms containing square-roots in the solutions to the field equations, and the requirement that solutions are real. From a more physical point of view, the singularity is related to repulsive effects originating from the Gauss-Bonnet term. Secondly, a minimum mass solution might exist in these kinds of theories, where the location of the finite radius singularity and of the event horizon overlap. In the above example, this critical solution possesses a non-vanishing Hawking temperature, presumably leading to the formation of a naked singularity.

III. EINSTEIN-SCALAR-GAUSS-BONNET GRAVITY: FIELD EQUATIONS AND THE SHAPE OF ξ (φ)
We now consider the more standard framework of the action of Eq. (2). Varying with respect to the metric tensor we obtain the Einstein equations where is the double-dual Riemann tensor (the square brackets denote anti-symmetrization). The scalar field equation is where the dot denotes differentiation with respect to the scalar field φ.
Starting with the scalar field equation (17), we review how different shapes of the coupling function ξ (φ) allow different phenomenologies. First we note that classical vacuum GR solutions require that φ = 0, and that these solutions only exists for couplings that obey the conditionξ(0) = 0. On the other hand solutions of models whose couplings obeyξ(0) = 0 necessarily differ from those of GR and possess a non-trivial scalar field. Common examples of couplings obeyingξ(0) = 0 are and The first, hereby dubbed the dilatonic (or exponential) coupling, is motivated from string theory, as it is the coupling that appears in the 4D low-energy limit of heterotic string theory [18][19][20][29][30][31][32][33]. The second -the linear couplingcan be considered as a linearization of the first around φ = 0, and additionally possesses a shift-symmetry in the scalar field [3][4][5]. For the dilatonic coupling, we focus on the γ = 1 case. Note, however, that as discussed in [55,57,58] the properties of black hole solutions might differ for other values of γ. In fact, for sufficiently negative values (γ −1), the small mass limit and the small size (radius) limit might not coincide, with small mass limit solutions being regular on and outside the horizon. For all the solutions discussed on this work, however, the minimum size and mass limits coincide, and therefore we will use these terms interchangeably. We chose to focus on the γ = 1 case because values γ −1 were explored in Refs. [55,57,58], and in our numerical explorations, all other values of γ to lead to similar behaviour as the one we observe for γ = 1. Therefore, we choose γ = 1 as a fiducial value. As discussed, if the coupling function satisfiesξ(0) = 0, then φ = 0 solves the field equations and the GR solutions are solutions to the model. If however,ξ(0) > 0 then the GR solutions are subject to tachyonic instabilities in the large curvature regime. This can be seen by linearizing the scalar field equation around GR solutions, for example about the Schwarzschild solution. For a perturbation δφ (for whichξ(0 + δφ) ≈ξ(0) +ξ(0)δφ) one finds that where G GR = 48M 2 /r 6 = 12r 2 H /r 6 > 0 (with M being the ADM mass of the black hole and r H the event horizon radius in Schwarzschild-like coordinates). Thus for α > 0, ifξ(0) > 0 the Schwarzschild black hole may develop an instability (as the effective mass gets negative, µ 2 ef f < 0) 4 . In this tachyonic regime, it has been shown that another class of solutions with a non-trivial scalar-field profile coexists with the GR solutions and are dynamically preferred, triggering spontaneous scalarization. In order to explore spontaneous scalarization we assume couplings of these type to obey the conditions The first condition can be imposed as the theory is invariant under ξ(φ) → ξ(φ) + cte, the second condition arises by requiring the existence of GR solutions and third condition can be imposed without loss of generality while maintaning a tachyonic instability in the large curvature regime. An example of such a coupling, commonly used in the literature [6,11,14,47,[68][69][70][71], and that we shall study here is the "quadratic exponential" coupling where β is a constant. Note that this choice is by no means unique, and a simple quadratic coupling ξ(φ) = φ 2 /2 (which is a particular case of Eq. (22) in the limit of vanishing β) would suffice to explore spontaneous scalarization per se. However, black hole solutions in models with a simple quadratic coupling are unstable, contrarily to those with a quadratic exponential coupling [9], and phenomenologies might differ. Note that the β = 3 case was studied extensively e.g. in Refs. [6,47]. We find that for any coupling obeying the conditions (21) the instability of a Schwarzschild black hole exists for (see Appendix A for a detailed discussion) Finally, we would like to point that another type of scalarization is possible, being induced by spin. For the Kerr metric, while for dimensionless spins χ ≤ 0.5 the Gauss-Bonnet term is positive definite, this is no longer true when higher spins are considered [72,73]. Therefore, if along with the other conditions in Eq. (21)ξ(0) is negative (instead of positive), fast-spinning Kerr black holes might be subject to a tachyonic instability [12,14,15]. A coupling compatible with this type of scalarization would be the coupling of Eq. (22) with a reversed overall sign. We will briefly discuss spin-induced scalarized solutions in Sec. VI.

IV. STATIC BLACK HOLE SOLUTIONS, AND THEIR SMALL MASS LIMIT
In the existing literature, it has been identified that static EsGB black holes exhibit similar behaviour to that described for the analytic case of Section II. Namely that in these situations there is also a minimum mass solution, beyond which solutions can no longer be described by black holes [4,18]. In the shift-symmetric case, it was further noticed that an inner singularity and the horizon overlap in this limit [4]. In this section we explore this small mass limit of EsGB black holes for a generic coupling function, discussing the domains of existence of solutions. To explore the small mass limit of static black holes in EsGB models we employ the static and spherically symmetric line element of Eq. (6), for which the field equations are presented in Appendix B. As already noted above no analytic solutions are known, and so numerical methods must be used.
Nonetheless, near the boundaries of our domain, analytic methods can be employed. We assume that a static black hole solution of the model allows the asymptotic expansion near the event horizon (hereby denoted by r H ) then, the near-horizon field equations to zeroth order in become which can be solved for f 1 and φ 1 in terms of r H and φ 0 , implying the following relations From these expressions one immediately finds the remarkable result that the horizon radius has a minimum size beyond which the solution can no longer be described by a black hole 5 . This follows by imposing the reality condition that f 1 and φ 1 must be real, leaving us with the condition or equivalently,ξ The above condition defines a region in the (r H / √ α,ξ(φ 0 )) plane within which BH solutions with a regular (real) scalar field configuration can exist. This is illustrated in Fig. 2, where we have also highlighted the spontaneous scalarization condition (23). Note, however, that the minimum horizon radius r min H depends onξ(φ 0 ). It is, therefore, possible that models where the minimum mass solution as given by Eq. (27) is never reached, depending on the behavior of φ 0 and of the coupling function. Indeed, as mentioned in the introduction, in Ref. [6] the authors follow the fundamental scalarized branch all the way to vanishing masses.
To further explore the small mass limit of EsGB black holes, consider again the field equations given in Appendix B, where we note that a closed-form expression for δ (r) can be obtained by a simple algebraic manipulation of the (r, r) equation. Substituting the value of δ onto the other field equations we can further rewrite the whole system of field equations in matrix form 6 where M and b are a 2 × 2 matrix and a 2 × 1 column vector respectively, whose components are given in Appendix B, and x = [f (r) φ (r)] T . Given the above system of equations, and appropriate initial conditions at some point r = r 0 , the existence theorem asserts that a solution to the Cauchy problem to extend our solution to a neighbouring point r 1 will exist if (see e.g. [53], which also studies the singularity structure of dilatonic black holes) Therefore if a point r = r * exists such that the determinant of M vanishes, the system of field equations will be illposed at that point. The existence of such a point would indicate the presence of a coordinate or physical singularity, whose nature would have to be studied by other means. From a numerical point of view, any standard strategy used for numerical integration will stop before r * . The determinant of M in Eq. (29) is given explicitly by where the expressions for A, B, C, D are again given in Appendix B. A simple examination of the determinant (31) reveals the existence of two zeros at the locations where f (r) = 0, and where The first, is related to the location of the event horizon (r * = r H ) and the correspondent singularity is a coordinate one. The second case is more intricate, and is related to a curvature (physical) singularity, as we will see. As a toy model to help us understand this singularity, consider Eq. (32) for a Schwarzschild background. This can be seen as the zeroth order solution in an expansion in α of the dilatonic and linear EsGB models. The solution is Similarly to the analytical example in Eq. (10), we observe that, at least for very small couplings, there is a singularity obeying (approximately) a proportionality relation r s ∝ (M α) 1/3 . Let us now explore the behavior of Eq. (32) near the event horizon of a EsGB black hole. Using the same near-horizon expansion as before (Eq. (24)), we observe that where r min H was defined in Eq. (27). Thus we observe that Eq. (32) vanishes at the event horizon if r H = r min H , indicating the presence of a singularity (other than the typical coordinate one). Therefore, we conclude that in the limit r H → r min H , an overlap of the curvature singularity and the event horizon occurs. In the following sections we will explore the small mass limit of EsGB black holes in more detail, utilising non-linear numerical solutions, and taking their inner structure into account.

A. Physical Quantities of Interest
To numerically integrate the field equations in order to obtain the black hole solutions, we match the near-horizon expansion of Eq. (24) with the appropriate asymptotic behaviour in far field (r → ∞) limit: where M is the ADM mass and Q s the scalar charge of the solution. Horizon quantities of physical interest include the Hawking temperature T H , the horizon area A H , and the entropy S. For the line element (6) these are given by where h is the determinant of the induced metric on the horizon, and R (2) its Ricci scalar [74]. The horizon and asymptotic quantities are related by a Smarr-type relation given by [74,75] M + M s = 2T H S, where For the dilatonic coupling (18) These conditions can be used to estimate the accuracy of our numerical method. Once again, an interesting remark can be made about the Hawking temperature, that in the small mass/size limit gives and thus, as in the case of section II, evaporation will not halt in the small mass limit and the black hole will continue to lose its mass at a rate given by Eq. (15), posing a potential threat to cosmic censorship 7 .  Fig. 4. We observe that detM has two zeros, one at the event horizon (blue dot-dashed line) and another at the singularity rs (red dashed line).

B. Numerical Method
We now compute numerical solutions to the field equations. We use a Runge-Kutta-45 ordinary differential equation solver and implement a shooting method for the parameter φ 0 such that the asymptotic expansions are matched. In more detail, the near horizon expansion of Eq. (24) is used to set initial conditions for a numerical integration, with the only free parameter being φ 0 (once r H and α are fixed). The field equations are then integrated from the horizon outwards to large r, the result is compared with the asymptotic expansion at large r, φ 0 adjusted, and the procedure repeated until the results of the numerical integration match the asymptotic expansion. Finally, using the results for the shooting parameters, the field equations are integrated from the horizon inwards to probe the internal structure of the black hole. We monitor curvature scalars such as the Ricci and GB scalars throughout the domain of integration, along with the determinant presented in Eq. (31). To test the accuracy of the numerical solutions we use the relations (38) and (40). We remark that errors are on the order of 10 −8 .

C. Numerical Results
Using the numerical algorithm described, we have explored the linear (19), dilatonic (18) and quadratic exponential (22) couplings. For the latter, we explore several values of β.

Linear and Dilatonic Couplings
For the linear and dilatonic couplings, by monitoring the Ricci and GB scalars, a finite radius singularity was always found at a radius r = r s > 0 inside the horizon, whose value ultimately depends on the ratio between the horizon radius and the coupling α. The singularity is located where detM in Eq. (31) vanishes inside the event horizon, as observed in Fig. 3.
In Fig. 4, we plot the metric functions along with the locations of r s and r H for several values of r min H /r H . As the horizon radius approaches r min H (as one would expect to happen dynamically as Hawking evaporation proceeds), the horizon and the singularity overlap, and numerical solutions reveal divergences of the curvature invariants, derivatives of the metric functions and the scalar field. That the location of r s and r H overlap when r H → r min H is in agreement with our analytical exploration shown in Eq. (34), and is similar to the behaviour observed in the analytical example of Section II (c.f. Fig. 1 and Fig. 5, right).
The domains of existence were constructed for both couplings, and can be observed in Fig. 5 (left), where each point on the dashed and dot-dashed lines represent a numerical black hole solution. Note that the domains of existence end at the red line, where r H = r min H . Assuming that Hawking radiation gradually reduces the mass of a black hole (and hence r H ) for some fixed α, but that our numerical solutions instantaneously remain an accurate approximation, we see that the fate of all black holes for both these couplings is to follow the lines on Fig. 5 and to always to reach the red line.

Quadratic exponential (spontaneous scalarization) coupling
Consider now the coupling of Eq. (22). We will perform a similar analysis as before, for several values of β. Note that higher values of β suppress scalarization. The domains of existence for β = 1, 3, 6 can be observed in Fig. 6 on the top left. We observe that for the β = 1 case, the domain of existence of solutions is similar to the dilatonic and linear coupling cases, where the location of the horizon and singularity overlap as the black hole shrinks, terminating in a critical solution. However, for β = 3 and 6, the domain of existence is radically different from the previous cases, as solutions never reach the singular (red) line, allowing the black hole to shrink to r H ∼ 0. This is possible due to the dependence of r min H on the value of the derivative of the coupling function at the horizon (as in Eq. (27)), and this behavior was observed for values β > β crit ≈ 2.33125, (42) while for lower values the black line would intersect the red one, the domain of existence ending in a critical solution.
The behaviour of the inner finite radius singularity is rather curious for β > β crit . We find there exists a finite radius singularity with r s > 0 only until to a certain value of r H / √ α, below which there is no singularity other than at the origin. This is shown in Fig. 6 (top right)   (bottom left), where the location of r s is plotted as a function of r H . On the bottom right in Fig. 6 we plot the location of the jump for a range of β. We have performed simulations for very large β of O 10 2 , and observed that the finite radius singularity exists only in a narrow of the domain of existence. For example, for β = 100 the finite radius singularity exists only from r H / √ α ≈ 0.83 (as in Eq. (23)) down to r H / √ α ≈ 0.828, and the maximum value of r s /r H is about 0.5.
On Fig. 7 we observe a fiducial numerical black hole solution (along with other relevant quantities) for which there is no finite radius singularity. Note that detM is strictly positive inside the event horizon. From a physical point of view we note that, as seen by the profile of the radial pressure p r in Fig. 7 (right), repulsive effects are maximum near the turning point of detM. The determinant then gets further away from zero as the repulsive effects get gradually weaker further inside the horizon.
An intuitive view on Hawking evaporation for this coupling would be the following. Starting from a (sufficiently large) Schwarzschild black hole, Hawking radiation will gradually reduce the mass of the black hole (for some fixed α), until the condition of Eq. (23) is met. A tachyonic instability would then settle in, leading to dynamical scalarization of the Schwarzschild black hole. The new scalarized solution will itself evaporate, with the endpoint now depending on the value of β. If β < β crit , the picture would be not too different to that of the dilatonic and linear couplings explored in the previous sections and Refs. [57,58], where the formation of a naked singularity is expected. However, if β > β crit , the evaporation process is expected to be similar to that of a Schwarzschild black hole, going all the way down to scales where quantum effects are expected to be important on general grounds, and our theory breaks down.
Remarkably, the behaviour of the inner singularity for the quadratic-exponential coupling shows that the critical solution end point of the domain of existence is not a generic feature of gravitational theories with higher-order curvature terms.

V. UPPER BOUNDS ON THE COUPLING CONSTANT
If evaporation proceeds as expected the behaviour we have been illustrating indicates that the small mass regime of EsGB theories may constrain the allowed form of couplings through self consistency arguments. However, even if evaporation is not taken into account, the minimum allowed size of black holes can place constraints on the strength of the allowed coupling through observational constraints.
To do so it is important to take a different form of the action such that results are consistent across the literature, and thus imposed on equivalent definitions of the coupling constant. When discussing observational constraints, the action for EsGB theories is usually presented in the following form [79][80][81][82][83] The mapping of the action (2) to the above parametrization can be done as and choosing ξ(φ) accordingly such that it is compatible with the definition of F (ϕ). To be consistent with the literature, we will constraint the coupling constant α.
For each of coupling functions, we consider the singular static solution with r H = r min H . Each of these singular black holes will have an associated minimum mass M min . From our numerical black hole solutions we can extract the quantity Assume now that an observation was made, in which a black hole was measured to have mass M obs . To be consistent with the description of a black hole within the EsGB theory we impose that the observed mass is greater than the allowed minimum mass Therefore, reintroducing the physical constants, from the above relation we obtain the bound This last equation allows us to impose an upper bound on the coupling constant for each coupling function. From our numerical solutions we have extracted Considering the case of GW190814 [84], where a compact object with a mass of around M obs = 2.6M was observed, and assuming it is a black hole, our calculations using Eq. (46) give the upper bound √ α 0.82 km, for To the best of our knowledge, these constraints are the tightest constraints on α so far, with the previous strongest upper bound being √ α 1.18 km for the linear coupling [79]. Constraints on the coupling obtained using data from other events can be found in Table I.  We have not studied constraints on the quadratic-exponential coupling in Eq. (22) for several reasons. First, there is an important dependence on β, i.e., if β > β crit , then as discussed above, no minimum mass solutions exist and no upper bound can be imposed. Also, in the small β limit we know solutions to be unstable [46,47]. Secondly, there is no guarantee that the black hole in question is a scalarized black hole.
Note that our results should be taken with a pinch of salt given that EsGB models are not UV-complete. It is possible that the pathological behaviour in the small mass limit is cured by higher-order corrections to the theory, depending on the scale they at which they become relevant.

VI. SPINNING BLACK HOLE SOLUTIONS
So far we have only studied static black hole solutions. In this section we will go one step further and explore spinning black hole solutions in EsGB theories. For this we resort, once again, to numerical integration of the field equations given in Section (III). The numerical procedure we adopt follows closely that of Ref. [5], Section 4. Namely, we consider a stationary and axi-symmetric line element of the form where the functions F i , W and the scalar field φ depend only on ρ and θ, and ρ H is the location of the event horizon in this coordinate system. The boundary conditions are implemented as follows. Asymptotic flatness is guaranteed by imposing while at the horizon, with the introduction of a new radial coordinate x = ρ 2 − ρ 2 H , we have where Ω H is the horizon angular velocity. Axial symmetry and regularity impose Focusing on black holes with parity reflection symmetry, we consider only the range 0 ≤ θ ≤ π/2, and impose the previous boundary conditions on the equator θ = π/2 instead of θ = π. The absence of conical singularities further imposes that on the symmetry axis F 1 = F 2 .
The system of coupled PDEs resulting from the EsGB field equations is solved with the FIDISOL/CADSOL solver [88][89][90], which implements a finite difference method together with the root finding Newton-Raphson method. We have also independently developed a new code to verify the accuracy of solutions which utilizes pseudo-spectral methods (and which will be described in a forthcoming publication [91]). The ADM mass M and angular momentum J of the black hole solutions are read of the asymptotic decay of the metric functions For future convenience we define the dimensionless spin of the solutions 8 Using the numerical method described above, in order to assess how the existence of a small mass limit changes with spin, we have explored the domain of existence of EsGB black holes for the linear, dilatonic, and quadratic exponential couplings.
For the dilatonic and linear couplings case (F (ϕ) = e ϕ and F (ϕ) = ϕ in Eq. (43)), the domain of existence of solutions with dimensionless spins χ 0.96 is displayed in Fig. 8 (left). The domain of existence is bounded by a critical line. As we approached the critical line numerically, we observed a divergent behaviour on the Ricci and Gauss-Bonnet scalars e.g. on the equator on the horizon, and the code eventually crashed. We observe that higher spins result in higher values of the minimum allowed mass (for fixed coupling). This is not unexpected, as spin adds another repulsive effect to the system. In turn, this translates to tighter upper bounds on the allowed value of α, if spin is considered. Therefore, the constraints of Table I constitute legitimate upper bounds. From the point of view of cosmic censorship and Hawking evaporation, black holes in both the linear and dilatonic theories are expected to give rise to naked singularities as their endpoint, regardless of the initial spin of the solution, again raising questions about the consistency of these theories altogether.
We now consider the spontaneous scalarization coupling of Eq. (22) with the particular value β = 6. In the static case we recall that no minimum mass was observed. We find, however, that when spin is brought into account, the picture changes and critical solutions do appear to exist, where curvature scalars diverge e.g. on the equator on the horizon, as observed in Fig. 8 (right). It is, however, unclear if the existence of rotating critical solutions changes the self consistency of the theory from the point of view of cosmic censorship and Hawking evaporation as it is not obvious that these solutions are ever reached. Indeed, a possibility is that during evaporation angular momentum is emitted at a (much) larger rate than mass, such that a rotating black hole spins down to a non-rotating state (which has no critical configuration) before most of its mass has been given up [92]. The endpoint of evaporation for this coupling is therefore an open question and constitutes an avenue of further research.
Finally, we have obtained preliminary results for spin-induced scalarized black holes, exploring several values of β (for the coupling of Eq. (22) with an overall reversed sign). In all cases, critical solutions were reached, in agreement with the results of Refs. [14,15].

VII. DISCUSSION AND CONCLUSIONS
In this work we have explored the small mass limit of stationary (both static and spinning) black holes in theories containing Gauss-Bonnet terms in the action. Starting with an analytical example, we explored the small mass limit of black holes in the theory known as gravity with a generalized conformal scalar field [59], which contains a Gauss-Bonnet term, and where static closed-form black hole solutions are known. These black holes do possess a minimum mass solution, where an inner singularity and the event horizon overlap. The inner singularity is intimately connected with the reality condition (that solutions must be real), because of the existence of terms containing square-roots on the solution (as is typical in Gauss-Bonnet theories). From a more physical point of view, the singularity is related to repulsive effects originating from the presence of the Gauss-Bonnet term in the theory.
Later, working with a more standard framework for EsGB theories, using numerical solutions of the field equations, a similar behavior was observed for the dilatonic (18) and linear couplings (19). A curious case concerns the quadraticexponential coupling (22), where for sufficiently high values of the constant β > β crit (defined in Eq. (42)), no static minimum mass solution was observed, thus showing that the existence of a critical singular black hole is not a generic prediction of theories containing Gauss-Bonnet terms. Then, from the point of view of cosmic censorship, this quadratic-exponential model might be viewed as more realistic option. The singularity structure for these models with β > β crit is rather different from that of the dilatonic and linear and merits a deeper study. Once spin is considered, critical solutions do exist, but it is unclear if these are ever reached from Hawking evaporation. Also, for the coupling of Eq. (22), scalarized black hole solutions exist only for curvatures above a certain threshold, rendering it particularly interesting.
Finally, we used the results concerning the minimum mass solutions into account to impose the tightest upper bounds to date on the coupling constant from observations, for both the dilatonic ( √ α (0.78 ± 0.03) km) and linear ( √ α (0.82 ± 0.03) km) theory, with the previous tightest upper bound being √ α 1.18 km [79]. Spin effects were found to only strengthen the previous upper limits.
where r H = 2M . Taking into account that the background geometry is static and spherically symmetric, the scalar field perturbation can be separated in the following way δφ = u(r) r e −iωt Y ,m (θ, ϕ), where Y ,m (θ, ϕ) are the spherical harmonics. The resulting equation for the radial part takes a Schrodinger-like form ( = 0) where dr * = dr/f (r) and A sufficient condition for the existence of an unstable mode (bound state as in quantum mechanics) is The above condition gives r H / √ α < 3/5 ≈ 0.774597 (or equivalently, M/ √ α 0.387298). Therefore Schwarzschild BHs with horizon radius obeying the previous condition, should be unstable in this framework. This, in turn, can be translated into a curvature condition: when the Gauss-Bonnet curvature at the horizon obeys the Schwarzschild BH should be unstable. This is only a sufficient condition for instability, but bifurcation of solutions actually occurs for slightly smaller curvature. To find the onset of instability we solve numerically the Schrödinger-like equation such that ω 2 = 0 (when ω 2 < 0 the tachyonic instability settles in). We find that the onset of instability occurs approximately at r H / √ α ≈ 0.83. This condition imposes a boundary on the domain of existence of (spontaneously) scalarized solutions. This result is valid for any coupling obeying the conditions of Eq. (21).

Appendix B: EsGB Field Equations for a Static and Spherically Symmetric Background
For a static and spherically symmetric background (21), the field equations take the form and the scalar field equation is where the primes denote a derivative with respect to r.