Black hole minimum size and scalar charge in shift-symmetric theories

It is known that, for shift-symmetric scalars, only a linear coupling with the Gauss-Bonnet invariant can introduce black hole hair. Such hairy black holes have a minimum mass, determined by the coupling of this interaction, and a scalar charge that is uniquely determined by their mass and spin for a fixed value of that coupling. Here we explore how additional shift-symmetric interactions affect the structure of the black hole, the value of the minimum mass, and the scalar charge.


Introduction
In recent years we have been presented with an unparalleled view of the strong field gravitational regime. The increasing number of gravitational-wave observations from the LIGO-VIRGO-KAGRA collaborations and the black hole images from the EHT collaboration have given us direct access to high-energy gravitational processes, including coalescing compact objects [1,2,3,4,5], accretion disks and shadows of supermassive black holes. The direct view of the highly dynamical, strong-field regime that current and future observations offer will give us an unprecedented chance to test the nature of gravity and search for new fundamental fields [6,7,8,9,10].
Indeed, black holes in general relativity (GR) are fully characterised by only their mass and spin and by the Kerr metric [11,12]. They can also carry an electric charge in principle, but this is expected to be entirely negligible for astrophysical black holes. Standard Model fields are not expected to endow black holes with any additional parameters [13,14], known as "hair", and hence any deviation from this picture in observations could imply the discovery of the new fundamental field or a breakdown of GR. Scalar fields have received particular attention in this context. There exist however several no-hair theorems which dictate that scalars cannot leave an imprint on quiescent black holes [15,16]. These theorems cover self-interacting scalars that potentially couple to the Ricci scalar and/or have a noncanonical kinetic term (scalartensor theories) [17,18]. However, more general nonminimal couplings, such as couplings to higher-order curvature invariants, are known to evade them, e.g. [19,20,21].
The case of shift-symmetric scalars is of particular interest when it comes to black hole hair, both observationally and theoretically. Strong gravity observations probe length (or curvature) scales of kilometres. Massive scalar profiles around compact objects are expected to decay exponentially and the characteristic scale for this decay is set by the inverse of the mass. Hence, these observations are expected to probe only ultralight or massless scalars. Shift symmetry, i.e. invariance under ϕ → ϕ+constant, is the symmetry that protects scalars from acquiring a mass. Consequently, strong gravity observations effectively probe scalars that exhibit either this symmetry or very small violations of it. Interestingly though, shift symmetry implies that the equation of motion of the scalar can be written as a conservation of a current, and this very property was used in [22] to prove a powerful no-hair theorem. It was subsequently shown in [23] that a linear coupling between the scalar, ϕ, and the Gauss-Bonnet invariant, G, evades this theorem and that this coupling term is unique in this regard. It is worth noting that the linear coupling term ϕ G would arise in a small coupling or small ϕ expansion of the exponential coupling e ϕ G that had already been known to introduce black hole hair [19,20,21].
Indeed, hairy black hole solutions in the ϕ G and in the e ϕ G cases share two key properties [24]. First, the scalar charge is not an independent parameter, but it is instead fixed with respect to the black hole mass and spin by a regularity condition on the horizon. Second, for any fixed value of the coupling constant that controls these terms (i.e. any given theory within the class), black holes have a minimum mass, controlled by the value of that coupling constant. The scope of this paper is to focus on the linear coupling case and determine how these two properties are affected by the presence of additional shift-symmetric (derivative) interactions in the action.
It is worth emphasising that having a minimum mass for black holes in the ϕ G model leads to a strong constraint on the coupling constant of this term, coming from the lightest black hole observed [25,26]. Most other observations are sensitive to the scalar charge but this gets converted to a constrain on the same coupling constant using the relation that fixes the scalar charge in terms of the black hole mass (and spin) e.g. [27,28,29,30]. In the case of extreme mass ratio inspirals (EMRIs), a scaling of the scalar charge with respect to the black hole mass inspired by the ϕ G theory proved to be a crucial ingredient to drastically simplify the modelling for a vastly broader class of nonminimally coupled scalar [29]. Note that, although there are good reasons to believe that such terms can remain subdominant when modelling binary dynamics and gravitational wave radiation [31,29], they could still have a crucial effect on the properties of the sources, including their quasi-normal ringing [32], and their dependence on the coupling constants of the theory [33]. This can affect how observations get translated to bounds for these coupling constants. Considering also that Effective field theory (EFT) strongly suggests that additional terms should be present, it is rather pertinent to understand their contributions to relations between mass and charge and to check whether they affect the minimum size of black holes.
The first step in this direction has already been made in [33]. It was shown there that, for the most general shift symmetric action that leads to second-order equations upon variation (shift-symmetric Horndeski theory) and respects local Lorentz symmetry, the scalar charge Q for black holes is given by 4πQ = α H n a G a , where H denotes the Killing horizon, n a its normal, G a in implicitly defined as G = ∇ a G a , and α is the coupling constant associated with the ϕG term. As expected from the discussion above, the charge vanishes if the ϕG term is absent. The value of G a does however depend on any additional couplings, with corrections with respect to the value G a would have in GR suppressed by the mass scales that correspond to these couplings. If one assumes continuity as these couplings are driven to zero and that their characteristic energy scales are similar to that of α (i.e. no hierarchy of scales), then one expects corrections to the charge and its scaling with the mass of the black hole to be subdominant. How much so is a matter of further exploration. Moreover, this expression for the charge does not give any information on the minimum size of black holes.
The structure of the rest of the paper is as follows: in Sec. 2, we will review the necessary theoretical background and present the class of theories we will consider. In Sec. 3, we will lay out the problem in a static, spherically symmetric setup and also consider the behaviour of a scalar field in a Schwarzschild background (decoupling) as a warm-up. In Sec. 4 we will derive hairy black hole solutions working perturbatively in the coupling α and analyse their properties, while in Sec. 5, we will present numerical results and a comparison with the perturbative ones. Finally, in Sec. 6, we present our conclusions.

Shift-symmetric Horndeski gravity and black hole hair
Horndeski's theory is the most general four-dimensional diffeomorphism-invariant theory involving a metric tensor and a scalar field that leads to second-order field equations upon variation [34,35]. We will restrict ourselves to shift symmetric theories. The shift-symmetric Horndeski action is then [24] where each sub-Lagrangian L i is given by is the Ricci scalar, and G µν is the Einstein tensor. We have also defined k = 8πG/c 4 with S M being the matter action. Matter is assumed to only couple minimally to the metric, that is, we are working in the so-called Jordan frame. Shift-symmetry implies that the field equation for the scalar can be written as a conservation of a current, The current is given by As discussed in the Introduction, this current conservation equation was used to prove a no-hair theorem in [22]. It was first shown that, for vacuum, static, sphericallysymmetric, asymptotically flat black holes, the only non-vanishing component of the current is the radial one J r . It was then argued that J r must vanish at the horizon, otherwise (J r ) 2 /g rr would diverge there. Current conservation then implied that J r must be zero everywhere. Finally, it was argued that J r = 0 everywhere implies that the scalar field must be constant everywhere. This no-hair theorem generalized to slowly rotating black hole straightforwardly [24]. It was however shown in [23] that a linear coupling between ϕ and the Gauss-Bonnet invariant, G = R µνρσ R µνρσ − 4R µν R µν + R 2 , circumvents this no-hair theorem. Indeed, G is a total divergence and so ϕG respects shift symmetry. However, consider the action The corresponding scalar equation is It can be written as a conservation of a current in the form ∇ µ (∇ µ ϕ + αG µ ) = 0, exploiting the fact that G is a total divergence. Although it is not obvious, the ϕG is part of action (1) [36]. However, it does not admit any constant ϕ solutions unless G = 0, which is not the case for black holes and hence they will have to have hair. The apparent contradiction with the no-hair theorem of [22] is resolved by the final argument of [22] -that vanishing current implies constant ϕ -relies on the assumption that every single term in the current depends on the gradient of ϕ. The contribution of the linear coupling to G clearly violates this assumption and is indeed unique in this respect [23]. Interestingly, hairy black hole solutions in this theory violate another assumption of the theory of [22]: (J r ) 2 /g rr diverges on the horizon [37]. It was subsequently shown in [38] that this quantity is not an invariant when it received a contribution from a linear coupling with G and hence there is no reason to impose that it is finite in this case.
In [33] theories described by action (1) where classified as follows: Class-3: All the rest.
Class-1 theories are defined by having ϕ = 0 as a solution for any general background; hence, they admit all possible GR solutions. Class-2 theories allow for ϕ = 0 to be realized only for flat spacetimes. The third class is defined as the complement of the other two. Therefore, class-3 theories admit a non-trivial scalar configuration in flat spacetime as a solution, or flat spacetime is not a solution, and hence they violate Local Lorentz symmetry.
At first sight, it appears that classes 1 and 2 are unrelated. On the contrary, it was shown in [33] that a class-2 Lagrangian can always be expressed as a class-1 Lagrangian plus a contribution from the Gauss-Bonnet invariant, namely: Consequently, all shift-symmetric non-Lorentz violating Horndeski theories admit all GR solutions, provided that a linear coupling between the scalar and the Gauss-Bonnet invariant is not present. Using this result, it was then shown that the scalar charge Q of a stationary black hole in any theory in classes 1 and 2 is given by where H denotes the Killing horizon and n a its normal, as already mentioned in the Introduction.

Our model
As already mentioned in the Introduction, the hairy black holes of the theory in action (7) have two key properties: their scalar charge Q is not an independent parameter, but it is instead determined by their mass (and spin), in accordance with eq. (13), and they have a minimum mass. In the next Section, we will see in detail how these properties relate to regularity conditions for static, spherically symmetric black holes. Our broader goal is to understand how adding additional shift symmetric terms to action (7) would affect these properties.
To make the calculations more tractable, we will not consider action (1). We will instead restrict ourselves to the following theory This action can be obtained from action (1) by selecting ‡ In units where G = c = 1 the scalar field is dimensionless while α, γ, σ, κ have dimensions of length squared. ‡ Up to a total derivative the term XR

Spherically symmetric setup
In this section, we consider a static and spherically symmetric background, described by the following metric, while the scalar field depends only on the radial coordinate, ϕ = ϕ(r).

Shift-symmetric current
The only non-vanishing component of the current J µ is J r , given by where a prime denotes a derivative w.r.t. the radial coordinate r. As discussed earlier, and according to the classification of [33], the current can be separated into a part for which every term contains ϕ ′ and a contribution by the coupling with G, as follows The conservation of the current can be straightforwardly integrated: Using (17), one can then determine ϕ ′ and then integrate once more to obtain ϕ(r). Due to shift symmetry, c in (19), is the only meaningful integration constant. Hence, considering also the mass parameter of the black hole, one would have a two-parameter family of solutions. The scalar charge would then be independent. However, ϕ ′ evaluated on the horizon of a black hole, r = r h , denoted as ϕ ′ h , generically diverges. If we assume that the scalar is regular on the horizon, and hence ϕ ′ h is finite, and we take into account that at r = r h we have A(r h ), B(r h ) → 0, then (17) implies thatJ r (r h ) = 0. Evaluating (18) on the horizon then fixes the value of c. As a consequence, the scalar charge ceases to be an independent parameter. After substituting c back in (18), and solving with respect toJ r we find: This is the equation we will be using in the following subsection.

Decoupling limit
As a warm-up, we consider the scalar in a fixed Schwarzschild background. We look for solutions that are regular at the horizon and approach a constant value asymptotically. For simplicity, we fix that value to zero, as the value of the constant is irrelevant due to shift symmetry. The γ term is not expected to contribute anything to the decoupled equations since it multiplies the Einstein tensor which vanishes at decoupling. The scalar equation on a GR-Schwarzschild background is: First, let us note that if only the α-term is present, one can find an analytic solution for the scalar field [23]: We note that, in this case, no specific restriction on the choices of α is suggested. We then consider the σ-term in addition to the α one, and we solve for the derivative of the scalar field: The quantity under the square root needs to be non-zero and positive, therefore, an existence condition emerges. It is straightforward to see that The inequality condition appearing in (23) imposes an upper bound on the product ασ. This, in turn, yields an upper bound on α when σ > 0, and a lower one when σ < 0. It is also possible to employ a near-horizon expansion, i.e. r = r h + ϵ, for the two cases discussed above. This yields When α is positive, from the near-horizon expressions we can deduce that for σ > 0 (σ < 0) the scalar field fall-off is larger (smaller) than in the σ = 0, α ̸ = 0 case, while in the limit σ → −∞ we retrieve the trivial solution ϕ ασ = 0 for the near-horizon expansion. When α is negative the aforementioned properties are reversed. The case where κ ̸ = 0 is more subtle. By examining the ϕ ′3 coefficient in (21), we see that when the κ-term is present, the derivative of the scalar at the horizon does not depend on κ, therefore, we deduce that κ does not enter an existence condition analogous to (23). If one attempts to solve the equation in the region [r h + ϵ, ∞), for ϵ ≪ 1, it turns out that a regular solution can be found ∀ κ ∈ R. However, not all of those solutions have the desired asymptotic behaviour, and for large positive values of κ, ϕ ′ ακ (∞) ̸ = 0.

Existence conditions
In the previous subsection, we saw how the existence conditions for the scalar equation are affected by the extra terms in the decoupling limit. Here, we derive the existence conditions for black holes beyond decoupling, for the full system of equations. To do so, we assume the existence of a horizon located at r = r h , so that A h ± → 0 ± , where the + sign corresponds to approaching the horizon from the outside, while the − to approaching it from the inside. By employing near-horizon expansions we obtain the following expression for the second derivative of the scalar at the horizon: We note that in order to get a black hole solution with a nontrivial scalar field it is required that (4αϕ ′ + r) ̸ = 0 [24] (cf. [20,39]). Since (A ′ /A) h diverges, in order for ϕ to be regular at the horizon, it is required that It is possible to derive from equations (26), and (27) two conditions: II: Condition I comes from the requirement that the quantity under the square root in equation (27) needs to be positive. Condition II comes from requiring that the denominator of the fraction in the right-hand side of equation (26) does not vanish on the horizon (where ϕ ′ is given by equation (27)). Although not obvious, these two conditions also guarantee that the denominator of the fraction on the right-hand side of equation (27) does not vanish.
It is worth pointing out that when we consider only the GB term, conditions I and II reduce to the existence condition appearing in [24]. Thus, we see that in the non-perturbative approach σ but also γ enter the existence conditions. In particular, in the case γ = 0, the parameter α has both an upper and a lower bound for either sign of σ. This will become more clear in sec. 5 where particular choices of the couplings are examined.

Perturbative treatment
First, we will employ a perturbative approach with respect to the coupling constant α, which is associated with the term that sources the hair. To do that we define the dimensionless parameterα ≡ α/r 2 h ≪ 1, where the horizon radius r h is the length scale we associate with our solution. In a similar manner we can defineγ = γ/r 2 h ,σ = σ/r 2 h andκ = κ/r 2 h . For nonzero, small values ofα we expect to acquire perturbative deformations to the Schwarzschild solution. Those are expressed through the following expansions for the metric elements Forα = 0, we retrieve GR minimally coupled to a scalar field, which for the spherically symmetric static configurations yields the Schwarzschild solution. These expansions are substituted in the equations of motion, which are then solved order by order for the unknown coefficients {A n , B n , ϕ n }. We work out the calculations up to the fifth order in the perturbative parameter O(α 5 ). The solutions become very lengthy beyond 2 ndorder and, therefore, are omitted. However, the expressions for the scalar charge and the mass of the black hole can be written in the following compact form: where the coefficients Q n and M n can be found in the appendix. Notice that we have to expand to the 5 th order inα before we see κ contributing to the scalar charge.
The perturbative treatment will break down at some radius. To trace when that happens we simultaneously scan the following expressions for perturbative inconsistencies. Note that the quantitiesĀ n ,B n appearing in equations (35) and (36) differ from A n , B n appearing in equations (30) and (31). If at some radius r np terms of different orders ofα become comparable in size, the perturbative treatment can no longer be trusted. The coefficients ϕ n , G n are given in Appendix B. From (35)- (38) we see that even at second-order inα, terms involving γ appear. We note that in the case whereγ,σ,κ = 0 it was shown in [24] that loss of perturbativity occurred at roughly the same radius at which the non-perturbative solutions exhibited a finite area singularity. We will return to this issue in the next section. In the top-left panel of Fig. 1 we present the radius r np denoting the point where the perturbative analysis breaks down, forγ = 0, ±0.5, ±1, ±5 andσ = 0. In the top-right panel of Fig. 1 we present theQ-M solution-existence curves for the same choices ofγ. The quantitiesQ andM are defined as: In the bottom panels, we present the analogous results for the caseσ = 0, ±0.5, ±1, ±5 andγ = 0. The horizontal axis in these plots corresponds to the normalized mass with respect to the GB coupling,M = M/α 1/2 , with M = 0.5. When r np exceeds r h part of the exterior cannot be described by the perturbative solution. One can read the corresponding mass from Fig. 1. From both top panels, we deduce that this mass increases/decreases for positive/negative values ofγ. We also see that forM ⪆ 2.5 all curves start merging, as the γ-term becomes significantly subdominant with respect to the GB one. One other interesting property we notice occurs for the value ofγ = 0.5 and it pertains to more than one solution existing for the same mass which can be seen in the right panel of Fig. 1. The radius of the singularity does not display similar behaviour and different mass black holes have different singularity radii, which is shown in the left panel of Fig. 1. These solutions are different from one another however as they describe black holes with different scalar charges.
From the bottom panels, we notice a similar trend regarding the effects ofσ and the mass for which perturbativity is lost already at the exterior. ForM ⪆ 3.5 all solutions in theM -Q plots begin to merge as the σ-term becomes subdominant. One of the main differences with respect to the top-panel plots, however, has to do with the maximum scalar charge. In theσ = 0,γ ̸ = 0 scenario we manage to get solutions with substantially larger scalar charges in comparison to theσ =γ = 0 (blue line). In theγ = 0,σ ̸ = 0 case this does not happen. Let us also point out that the values for r np depicted in the bottom left panel, are the same for positive and negative values of σ. To understand why this occurs it is helpful to see the way σ enters the perturbative expansions, which are presented in Appendix B. Specifically, σ appears at second order in the expansion of the scalar field as a multiplicative constant, which explains why changing its sign yields the same solution for r np .
It should be noted, however, that these conclusions have to be drawn with care, as they correspond to a region of the parameter space whereM < 5, orα > 0.01. This is a region that can in principle render the perturbative approach problematic in general and only a proper numerical analysis can either confirm or disprove the aforementioned effects.

Numerical results
We now move to solve the full system of equations numerically. This is a system of ordinary partial differential equations (ODEs) of the form {ϕ ′′ , A ′ , B ′ } = f (r, ϕ ′ , ϕ, A, B). We separate the analysis into two regions: the black hole exterior and the black hole interior. In both cases, the integration starts at the horizon. The theoretical parameter space consists of (γ, σ, κ, r h ), where r h is the black hole horizon radius. Since r h appears in the existence condition (27), the allowed values for the coupling parameters are expected to be affected if we variate r h . We can straightforwardly reduce the dimension of the parameter space by one if we normalize the coupling parameters with the horizon radius as we did in the previous section.
For a given theory defined by (γ,σ,κ) we allow the values ofα to scan the parameter space starting from smallα and gradually increasing until the existence conditions are saturated, We, therefore, need the set of values {ϕ ′ , ϕ, A, B} r h . Despite appearing to constitute "initial data", this set of values is not entirely free to choose. In practice, in order to apply the existence conditions (28)- (29) with reasonable numerical accuracy, we use a perturbative expansion near the horizon and we numerically solve the system of algebraic equations for the first few coefficients appearing in the expansions (up to order O(r − r h ) 2 ). This process reduces the number of the free initial conditions to two, namely the value of the scalar field at the horizon and that of the first-order coefficient of A. The latter one, however, is fixed by asymptotic flatness, leaving ϕ h as the only free-to-chose initial condition. The asymptotic value of the scalar field should be constant but otherwise unconstrained since our model is shift-symmetric. For simplicity, we choose ϕ h so that ϕ ∞ = 0. To achieve that we employ a shooting method while integrating outwards, demanding ϕ vanishing to a part in 10 4 . The remaining free parameter is the horizon radius r h . We then start the numerical integration outwards (inwards) from r = r h ±O(10 −5 ). In the exterior, we typically integrate up to r/r h ≈ 10 5 .
In the following subsections, we present plots corresponding to different cases of couplings. In each case we numerically calculate the scalar charge and the the Arnowitt-Deser-Misner (ADM) mass of the black hole using the following expressions: We were also able to verify the emergence of a finite-radius singularity, consistent with the existence conditions. While integrating from the horizon and inwards we noticed the following general trend: Starting from GR (α → 0) and gradually increasing the couplings, the geometric invariants diverge and the solutions become singular at some radius r s . The larger GR deviations become the more r s approaches r h . When one of the existence conditions is saturated the singularity radius approaches the horizon radius, i.e. r s → r h .  In these subsections we attempt to present the overall generic trend that the black hole properties follow, if one considers action (14). We discuss the charge, mass and scalar profile for a few examples corresponding to different scenarios, which motivate the more thorough analysis that follows in the next subsections.

Charge, mass and scalar profile
On the left panel of Fig. 2, we show theM -Q plot for different negative values of the coupling constantsγ andσ. The corresponding plots for positive couplings are not presented here, since -at these scales-they are overshadowed by theγ =σ = 0 curve, as will be explained later in more detail. In all positive-coupling cases, the minimum black hole mass is larger than the one corresponding toγ =σ = 0. Furthermore, non-zerõ κ curves are almost indistinguishable from theγ =σ = 0 one since as we sawκ does not enter the existence conditions. Consequently, the correspondingκ-plots are not presented here. We see that for largeM which corresponds to small GB couplings, the charge in all cases drops off to zero and GR is retrieved. This is of course associated with the fact that the GB term is the one sourcing the hair. In the smallM regime significant deviations are observed, which are explained in the following coupling-specific subsections. On the right panel of Fig. 2, we show the profile of the scalar field, properly normalized with the distance and the scalar charge. All curves exhibit a 1/r fall-off and asymptotically approach 1. For large radii, the scalar field profiles are indiscernible for different couplings. In the near horizon regime, however, there are apparent deviations in accordance with the non-trivial deviations shown in the left panel.
To make things easier for the reader, in what follows, we consider the GB coupling α in combination withγ,σ andκ separately. In this work, we consider α > 0 as this is consistent with most of the bibliography. However, it is worth pointing out, that action (14) is invariant under the simultaneous transformation α → −α, ϕ → −ϕ and σ → −σ and that in the case of σ = γ = κ = 0, the sign of ϕ is determined by the sign of α for solutions that are continuously connected to Schwarzschild as α → 0. In what follows we consider both positive and negative values for σ, γ, and κ, and hence our analysis should effectively cover the α < 0 case as well, at least for configurations that are continuously connected to Schwarzschild.

Theγ term
First we consider the caseσ =κ = 0. From (28) and (29) we find the conditions oñ γ necessary for regularity at the horizon. The existence conditions I-II are in general non-trivial and the easiest way to track them is to examine the corresponding region plot. In the left panel of Fig. 3 we see the aforementioned plot withγ being on the horizontal axis andα occupying the vertical one. The first obvious observation relates to the apparent asymmetry about the vertical axis. Therefore, we expect the sign ofγ to influence the black hole solutions and properties. In particular, for negative values ofγ the parameter space of allowed values forα increases, and so we expect negative values ofγ to allow for hairy solutions with smaller masses. On the other hand wheñ γ > 0 the parameter space ofα shrinks and we expect the black hole mass range to also decrease. Regarding the GB-couplingα, extending the plot to negative values ofα is trivial as, forσ = 0, the action (14) is invariant under the simultaneous transformatioñ α → −α and ϕ → −ϕ.
These are indeed verified in Fig. 4, where the emergence of a finite-radius singularity is demonstrated in the interior of the black hole. The left panel shows the singularity radius of the black hole mass forγ = {0, 0.5, 1, 5} while the right panel shows the corresponding results forγ = {0, −0.5, −1, −5}. The values are chosen to be of order ∼ 1 − 10 with respect toα max , whereα max corresponds to the largest allowed value for α satisfying the existence conditions. For the choices ofγ made, we present the results for the minimum hairy black hole mass in the following table: Minimum mass forσ =κ = 0,γ ̸ = 0,α > 0 For a negativeγ, we notice another interesting property of the solutions: at small masses, the apparent change in monotonicity in theM -Q andM -r s (see the inset) plots indicates that black holes with the same mass can correspond to different scalar charges and singularity radii. Therefore, one would expect that the black hole with the larger scalar charge, would shed some of it to reach a more favourable scalar configuration with a smaller charge. Finally, from Fig. 4 it is pointed out that in the larger mass regime, the sign ofγ becomes unimportant and the cases with opposite signs merge.

Theσ term
In the left panel of Fig. 5 we present the allowed and excluded regions of the parameter space according to the existence conditions in the case ofγ = 0, withσ being on the horizontal andα on the vertical axis. For negative values ofα the region plot we retrieve demonstrates an origin symmetry which was anticipated since the action (14) is invariant under the simultaneous transformationα → −α,σ → −σ, and ϕ → −ϕ. Forα > 0,σ < 0 the allowed values forα increase and therefore the mass range also increases, and hairy black holes with smaller masses are found. At the same time for α > 0,σ > 0, the parameter space ofα shrinks and the black hole mass range should also decrease. If we consideredα < 0 the above conclusions would be reversed. In Fig. 6 we display the singularity radius in this scenario and its dependence on the value ofσ. Verifying the above, positive(negative)σ leads to a larger(smaller) minimum black hole mass. In theσ ̸ = 0 scenario, the relation between the finite singularity radius and the normalized mass exhibits discontinuous behaviour, which is evident from the vertical jumps shown in Fig. 6. As already explained, we identify the singularity radius as the one for which a geometric invariant (e.g., the Gauss-Bonnet or equivalently the Kretschmann invariant) diverges. To explain the discontinuity let us imagine that we start from some largeM moving inwards towards smaller masses. At r = r s the GB invariant diverges and we identify r s as the singularity radius. There exists a second special point at r = r ′ s > r s where the metric functions and the scalar field appear to lose differentiability because they develop a kink. This appears to be because the second derivative become discontinuous. The differential equations however can still be integrated for r ′ s > r > r s . If we plotted r ′ s instead of r s , then the vertical jump would no longer be present and the lines would be continuous. In all cases, however, we chose to plot the singularity corresponding to the divergence of the geometric invariants. On the other hand, for positiveσ, we do not encounter any other "singularities" than the ones we plot, which correspond once again to the geometric invariants diverging.
Similar discontinuities have also been encountered in scalar Einstein-scalar-Gauss-Bonnet gravity with a quadratic exponential coupling [25]. In the zoomed-in part of the right panel of Fig. 6, similar behaviour to the negativeγ case is exhibited, where samemass black holes have different singularity radii. This can also be understood from thê M -Q plot, in the right panel of Fig. 5, where a turning point appears at small masses.

Theκ term
It is evident from equations (28) and (29) that κ does not enter the existence conditions. As a result, one might naively conclude that black hole solutions exist irrespective of the value that κ takes, given that the remaining parameters satisfy the existence conditions. Contrarily, that is not the observed behaviour. If κ is taken to be positive, then we cannot find solutions for all values of α that are allowed by the existence conditions; however, if κ is negative, then solutions could be found for all values of α allowed by the conditions. This behaviour is illustrated in Fig. 7, where for negative values of κ it is possible to saturate the existence condition and have solutions with a naked singularity, but for κ > 0, in general, that cannot be achieved. To better understand this trend, it is useful to rewrite the scalar equation for γ = σ = 0 as In practice, we see that whenκ > 1 not all values ofα, allowed by the existence conditions, yield black hole solutions. In order to give an explanation to this issue, we numerically examine the value of the quantity inside the square brackets in eq. (41),  namely h µν . Due to the symmetry of our problem only h rr will be examined. In Fig. 8, we plot h rr (and h r r ) for values ofκ spanning a few orders of magnitude, i.e. O(1) − O(10 2 ). We see that for O(κ) > O(1), the quantity h rr approaches zero at some intermediate radius, which seems to increase as we increase the value of κ. Beyond that point, the ODE system can no longer be integrated. This bares similarities with the behaviour of ϕ ′′ at the horizon, the regulation of which yielded the existence conditions. Thus, it appears that imposing regularity for the scalar field at the horizon results in divergences appearing elsewhere for large positive values ofκ.

Numerical solutions vs perturbative solutions
As mentioned earlier, it has already been demonstrated that in the case γ = σ = 0 loss of perturbativity is associated with the appearance of a finite-radius singularity in the black hole interior. Here we discuss the relation between the perturbative treatment breakdown radius r np and the finite-radius singularity r s in the general case whereγ,σ are nonzero. We present the comparative plots in Fig. 9. Verifying the results of [24], we see that the radius of the singularity in the black hole interior in the caseγ =σ = 0, is traced almost perfectly by the perturbative analysis. However, this is not the case, at least to the same level of success, when one considers the γ and σ contributions. From the left panel of Fig. 9, we see that whenγ ̸ = 0 the r np curve sits below the singularity radius r s . From the right panel, we notice that in theσ ̸ = 0 case on the other hand the r np curve sits between the disconnected branches of the numerical solutions. We can also compare the results regarding the scalar charge and mass of the black holes, by inspecting Figs. 1 and 3, 5. Specifically, for theγ ̸ = 0 case, by comparing Figs. 1 and 3 we deduce that the perturbative approach captures at least quantitatively the two main conclusions drawn by the numerical analysis: (i) Forγ > 0 the minimum black hole mass for eachγ increases and the mass-parameter space of solutions shrink. (ii) Forγ > 0 the minimum black hole mass for eachγ decreases, and for small masses it is possible to retrieve black holes of the same mass with different scalar charges. Similar conclusions were also drawn in the numerical analysisσ ̸ = 0 scenario (with the addition of the discontinuities), where the qualitative trends were also captured by the perturbative analysis.

Conclusions
We have studied hairy black holes in generalized scalar-tensor theories that exhibit a range of shift-symmetric derivative interactions, in addition to the linear coupling to the Gauss-Bonnet invariant that is known to introduce black hole hair. We found that, although these additional interactions cannot introduce hair themselves, they can significantly influence the behaviour of the scalar fields near the horizon of the black hole and hence affect the configuration in general, including the value of a scalar charge for a given mass.