Black holes in self-tuning cubic Horndeski cosmology

Observations of neutron star mergers in the late Universe have given signiﬁcant restrictions to the class of viable scalar-tensor theories. In this paper we construct black holes within the “self-tuning” class of this restricted set, whereby the bare cosmological constant is absorbed by the dynamics of the scalar, giving a lower eﬀective cosmological constant. We use analytic expansions at the singularity, black hole and cosmological horizon, and asymptotic region, coupled with numerical solutions, to ﬁnd well-behaved black holes that asymptote to the self-tuned de Sitter geometry. The geometry diﬀers from standard general relativity black holes near the horizon, and the scalar ﬁeld velocity provides a hair for the black holes.


I. INTRODUCTION
It is a remarkable success of the theory of General Relativity (GR) that it has so far withstood the onslaught of ever more precise astrophysical and cosmological data, and at present remains our best low-energy description of gravity. That is not to say, however, that GR is without its problems, particularly in providing a robust description of the dark sector of the universe. These outstanding issues have motivated research into modifications of GR, the prototypical example being scalar-tensor theory, in which one introduces an additional scalar degree of freedom into the gravitational sector. These are typically described by Horndeski theory [1][2][3][4], or its more recent generalisation, Degenerate Higher Order Scalar Tensor (DHOST) theory [5][6][7]. Such an approach has proven particularly popular in recent years for constructing novel descriptions of dark energy. The recent gravitational wave data [8][9][10][11][12][13] has dealt a strong blow to such models, proving that in the late Universe light and gravity propagate at identical speeds, and thus ruling out large regions of the parameter space of these theories. However, a simple subclass of Horndeski Lagrangians passes this test. Some models even provide an interesting self-tuning mechanism.
Self-tuning is the ability of a model to screen an arbitrary cosmological constant at the level of the action, and possibly to yield an independent, effective cosmological constant at the level of the solutions. Examples of self-tuning towards Minkowski spacetime were first provided in [14] and dubbed 'Fab Four' models; the mechanism was later generalized to non-zero effective cosmological constants (e.g. [15][16][17]). This is particularly interesting in the context of the large cosmological constant problem; as is the case for any parameter in a quantum field theory, the cosmological constant that we observe should be the sum of a bare contribution and of quantum corrections. However, the bare value of the cosmological constant and the associated quantum corrections must compensate at a level of at least one part in 10 55 [18] in order to yield the observed value. Self-tuning is thus extremely useful, as a way to trade this huge contribution for an effective cosmological constant that matches the observations.
On the other hand, with gravitational wave detectors and electromagnetic observations, we have entered an era that makes it possible to explore the strong field regime of gravity. It is thus more important than ever to understand the behaviour of scalar-tensor theories in this sector. Of particular interest is the study of astrophysical black holes, to which current observational and detector experiments are sensitive to. Indeed, the incoming data provides a new testing ground in which to probe the validity of such modifications to GR. Determining black hole solutions within scalar-tensor theories is, however, a complicated endeavour, due to the presence of non-minimal couplings between curvature terms and the additional scalar degree of freedom in the gravity sector. Nevertheless, the problem of embedding black hole solutions in self-tuned universes has already been tackled. Such solutions have been obtained through approximate [19], numerical [20] or even exact methods [21][22][23][24]. Star solutions were also obtained numerically in similar contexts [25][26][27][28]. Exact black hole solutions were only obtained in the framework of shift-symmetric theories, when the action is invariant under an overall shift of the scalar field across spacetime. In this paper, we also focus on the subclass of shift-symmetric theories; this allows us to maintain tractable calculations, while capturing the higher-derivative effects that are essential to self-tuning.
The structure of the paper is as follows. In sec. II, we introduce our model and the associated field equations. Then, in sec. III, we discuss its self-tuning properties. We present the only self-tuning model left in the class of Horndeski theories that pass the gravitational wave tests, together with an exact cosmological solution. Section IV describes how black holes can be embedded in such a universe, both using analytic expansions and numerical techniques. We conclude in sec. V.

II. A SHIFT-SYMMETRIC SELF-TUNING MODEL
The subset of Horndeski theory that passes the gravitational wave test can be parametrized as [8][9][10][11][12][13] where M Pl is the reduced Planck mass, Λ is the bare cosmological constant, and K, G 3 and G 4 are a priori arbitrary functions of the scalar field φ and its kinetic density X = − 1 2 ∇ µ φ∇ µ φ. Note that G 4 is taken to be a function of φ alone, in order for the theory to remain consistent with the gravitational wave constraints. The element that is essential to self-tuning is the X dependence in G 3 . Thus, for tractability of the calculations, we will consider the simpler shift-symmetric subset of (1): which results in the following action: The above action was actually first studied in [29], where it was dubbed Kinetic Gravity Braiding (KGB). It was later incorporated in the generic framework of covariant Galileons [2][3][4], which were in turn found to be an equivalent formulation of Horndeski's theory [1]. The associated metric and scalar field equations of motion are then given by: where is the conserved Noether current associated to the shift-symmetry, and the subscript X denotes a partial derivative with respect to X (e.g. G 3X = ∂G 3 /∂X).

III. COSMOLOGICAL SOLUTIONS TO THE FIELD EQUATIONS
Let us first study the cosmological solutions of the model, within which we aim to embed black holes. Any black hole solution should match this behaviour asymptotically. We start with the Friedmann Lemaître Robertson Walker (FLRW) geometry: where τ is the cosmic time, a(τ ) is the cosmological scale factor, and we adopt the metric signature (−+++) throughout. The scalar field is also assumed to be spatially homogeneous at this scale, i.e. φ = φ(τ ). Given this, the metric and scalar equations of motion, eqs. (4)-(5), are where a dot denotes a derivative with respect to cosmic time τ , H =ȧ/a is the Hubble parameter, and G 3 depends implicitly onφ through X =φ 2 /2. We now fix the geometry to be de Sitter, i.e.Ḣ = 0, H = H 0 . It immediately follows from eq. (8a) thatφ is also constant. The G 3 dependence then drops from eq. (8b), from which we can determine thaṫ φ is given byφ Once substituted in eq. (8a), one realises that H 0 will always depend on Λ, unlessφG 3X is itself insensitive to the value ofφ. Thus, to ensure self-tuning, i.e. that the cosmological properties of the metric do not depend on the bare cosmological constant, we set where µ is a priori an arbitrary parameter in the action (1), with mass dimension [µ] = 1. Interestingly, the particular choice of G 3 in eq. (10) corresponds to a constant diffusion coefficient in the imperfect fluid interpretation of the KGB model [30]. The on-shell field equations then impose that H 0 = µ, such that the value of the Hubble parameter is independent of the cosmological constant, i.e. the model admits self-tuning solutions. Having to introduce the current Hubble scale of µ ∼ 10 −33 eV directly at the level of the action can appear unnatural. However, it is counterbalanced by the advantage that H 0 is completely independent on Λ; thus, the tuning will survive phase transitions, that are otherwise extremely problematic in the context of the standard model [18]. We have established the cosmology to which we will match black hole solutions. We now need to consider how we can achieve such a matching. We will look for static and spherically symmetric black hole solutions. Thus, we will use the following system of coordinates: where dΩ 2 = dθ 2 + sin 2 θ dϕ 2 , and f (r) and g(r) are a priori arbitrary functions to be determined by the field equations 1 . The de Sitter metric, eq. (7) with a(τ ) = e H 0 τ , can actually be matched to its well-known static patch through the following invertible change of coordinates [19]: Under this change of coordinate, the metric (7) with a(τ ) = e H 0 τ is mapped to IV. STATIC BLACK HOLE SOLUTIONS

A. Ansatz and field equations
Having determined how a black hole solution within this shift-symmetric model should behave asymptotically, i.e. on cosmological scales, we shall now move on to determine the details of such solutions. For simplicity, we shall consider static and spherically symmetric metric configurations, as given by the ansatz (11). The scalar field, on the other hand, is taken to depend not only on the radial coordinate r, but also linearly on time. This emerges naturally from the combination of eqs. (9)-(12a), and was shown to be fully consistent in [31]. The stress-energy tensor associated with φ only involves derivatives of the scalar field; thus it is itself static and one can consistently require a static geometry 2 . Furthermore, although a naive counting of the field equations and free functions tells us that the system seems overconstrained, ref. [31] proved this is not the case; the (tr) component of the metric equation is always proportional to the radial component of the current, J r . This allows us to look for solutions with the following ansatz for φ: where q =φ is a free, non-zero constant for now, and χ(r) is an arbitrary function to be determined by the field equations. The precise form of the second term on the right hand side is chosen for convenience. Notably, the f (r) denominator allows us to absorb the diverging contribution of the r-dependent part of φ when approaching the black hole horizon. Let us now write down the independent field equations. We will work with the (tr) component of the metric field equations. Additionally, we use specific linear combinations of the (tt) and (tr) field equations on one hand, and (rr) and (tr) field equations on the other. Due to spherical symmetry, E θθ (g) = sin 2 θ E ϕϕ (g) , and E θθ (g) itself can be deduced from the previous three equations, due to the diffeomorphism invariance identity E (φ) ∇ ν φ + 2∇ µ E µν (g) = 0. The (tr) metric equation reads (up to an overall non-zero factor): with Then, the (rr)-(tr) and (tt)-(tr) combinations read respectively: Note that the G 3 dependence drops out from the (rr)-(tr) combination, eq. (17). It is not possible to solve this system fully analytically, and therefore, we will ultimately have to appeal to numerical techniques. The difficulty of finding exact black hole solutions for cubic Horndeski models was already noted in [20]. In general, it appears that no exact black hole solution is known for Horndeski models that do not possess the reflection symmetry φ → −φ, particularly the simplest cubic and quintic models. Indeed, typically the scalar field equation can be integrated once, but it then becomes a high-order algebraic equation for φ . Exact solutions are known only for models that possess both shift and reflection symmetry [21,22,32] (apart from standard scalar-tensor theories). In order to integrate the system numerically, it proves convenient to introduce a set of dimensionless quantities.
To this end, we note that the Lagrangian contains three dimensionful parameters: M Pl , µ and Λ, along with the scalar field velocity q, introduced in the ansatz for φ. Let us further introduce a length scale r 0 , which physically corresponds to the event horizon radius of the black hole, such that we can define a dimensionless radius x = r/r 0 . Then, given the parameters of the Lagrangian, the scalar field velocity q and the length scale r 0 , we can express the field equations using only the following three dimensionless constants: Equations (15), (17) and (18) can then be rewritten: These will be the equations that we integrate numerically.

B. Analytical approximations in the small & large r limits
Before solving eqs. (20), (21) and (22) numerically, we shall study the asymptotic behaviour of f , g and χ in the small and large r limits. Solving the system of equations near the origin r = 0, we find that where a, b and c are fixed by the field equations (15), (17) and (18) (we omit their exact expressions here in the interest of succinctness). We note that, as was found in [20], the g(r) component of the metric is finite at the origin, unlike in GR. It is possible that this is a generic feature of any cubic Horndeski model for which such black hole solutions exist.
Let us now consider the large r asymptotic behaviour of the solutions. We assume that, as r → ∞, the solution has an analytic expansion in powers of 1/r. We find that the asymptotic behaviour of the solution is given by: where q 0 =φ 0 is given by eq. (9). Note that, at this point, the velocity parameter q remains arbitrary, and may not necessarily coincide with q 0 , i.e., the velocity parameter corresponding to the cosmological solution. However, it is interesting to check whether the scalar field in the asymptotic solution (24) is homogeneous in the set of cosmological coordinates (τ, ρ). From eqs. (14) and (24), it follows that Similarly, using the coordinate transformation (12), φ can be expanded as in cosmological coordinates. Thus, we see that φ is asymptotically homogeneous only when the local velocity q matches the asymptotic cosmological value, as was already the case in [20]. Whether the solutions with q = q 0 are physically acceptable is arguable, since only the gradient of φ is physically relevant, and the latter always decays as 1/ρ. In our numerical analysis, still, we chose to impose that q = q 0 to ensure a safe cosmological behaviour. Finally, let us comment briefly on the fate of the scalar field at the (black hole and cosmological) horizons. It is well known [21] that in the system of coordinates (11), the rdependent part of φ diverges at both horizons. However, this is a coordinate effect, and this divergence can be taken care of by working in adapted Eddington-Finkelstein coordinates. Note that it is not in general possible to absorb the divergence for both horizons. However, recent works have shown that generalising to the rotating case can cure the divergence of φ on both horizons for this type of solution [33].

C. Numerical integration of the field equations
Having studied the asymptotics of the solutions, we shall now proceed to carry out a numerical integration of this system. Note that eqs. (20) and (21) are algebraic equations in g and χ, and can thus, in principle, be solved to express them in terms of f and df /dx. Upon substitution into eq. (22), we are then left with a second-order ordinary differential equation (ODE) for f . In practice, it turns out the system is more readily solved numerically by first solving eq. (21) to determine χ in terms of f , df /dx and g, and then solving the remaining two equations as a system of ODEs for f and g. This requires us to specify several boundary conditions in order to obtain the unique solution. By imposing that both f and g vanish at x = 1, we define the black hole horizon to be located at this point (corresponding to r = r 0 as mentioned earlier). Moreover, having imposed f | x=1 = 0, we can then expand eq. (20) around x = 1 to determine df dx | x=1 in terms of dg dx | x=1 (or vice versa). This leaves us with one remaining boundary condition that we are free to choose arbitrarily, although ultimately the choice is dictated by requiring that at large x, the solution has the desired asymptotic cosmological behaviour. Thus, for given parameters, we use the shooting method to pick the value of dg dx | x=1 such that f ∼ x→∞ g. We start the numerical integration slightly outside the black hole horizon. Then, the code breaks down when approaching both the black hole and cosmological horizons. It is however very easy to perform a series expansion around these points, and to restart the numerical integration on the other side of the horizon(s). Thus, we are able to compute the solution over the full range 0 < x < +∞. A typical solution is shown in Fig. 1 In this plot, the parameters are chosen to be β 1 = 100, β 2 = 1.76, β 3 = 1. This corresponds to a modest hierarchy of β 1 β 3 = 10 2 between the bare and effective cosmological constants. The ratio can be increased, but we kept these scales for clarity of the plots. and 3 respectively. The numerical solutions agrees with the analytic expansion presented in sec. IV B. As mentioned in sec. IV B, we stick to q = q 0 in order to ensure the good asymptotic behavior. In terms of dimensionless parameters, this translates as The metric and scalar field functions at large radius, together with the asymptotic expansion of eq. (24). Again, the parameters are the same as in Fig. 1, and the convention in the legend is the same as in Fig. 2.

V. DISCUSSION
We have considered the class of Lagrangians within Horndeski theory for which light and gravity travel with identical speeds, as motivated by observational data in the late Universe [8][9][10][11][12][13]. We restricted ourselves to shift-symmetric models, in order to keep the essential higher derivative effects while handling reasonably tractable equations. In this context, only one model exhibits a fully self-tuning behaviour, allowing an exact de Sitter solution fully independent of the bare cosmological constant. We discussed how black hole solutions can be matched with these cosmological asymptotics using an analytic expansion. We then explicitly built such solutions using a numerical shooting method. The numerical solutions agree with the analytic expansions, both at large radius and close to the central singularity. We notice that, at the centre of the black hole, one of the metric components remains finite, a behaviour that was already found in [20] for a similar model.
It is interesting to note that the Lagrangian (3)-(10) falls in the well-tempered category [34,35]. The well-tempered proposal is somehow similar to self-tuning, but it aims at treating dark energy separately. The motivation of the well-tempered model is to resolve a problem of the previous self-tuning "fab-four" model [14], whereby the model's ability to absorb a huge cosmological constant came at the price of it greatly affecting any matter/radiation content as well. A consistent cosmology requires, in particular, a radiation and a matter dominated era, so one should be careful to allow such an era in the cosmological history. In the well-tempered set-up, dark energy is on a special footing due to its particular equation of state P + ρ = 0. This allows for the scalar field to only self-tune for constant energy densities, such that radiation and matter dominated eras seem to be generically present. However, in this paper, we focused on exact de Sitter solutions for the asymptotics. In the well-tempered context, exact de Sitter solutions exist only for a certain relation of the type H µ √ Λ, which we obviously want to avoid. Thus, although the model we study is in the well-tempered class, we do not focus on a well-tempered branch but on a self-tuning one. Still, it seems that well-tempered models possess solutions that approach H µ √ Λ asymptotically [34,35]. An interesting extension of our work would be to study the welltempered branches. This however comes with the extra difficulty that the solution is only approximately de Sitter at the cosmological level. It cannot be matched to a static patch, and one would likely have to consider non-static metrics.
We also note, following ref. [24], that local gravity tests (such as Solar System tests) are likely to fail for our model if one wants to tune too large a value of the bare cosmological constant. This is true even when non-linear effects, usually responsible for a Vainshtein screening, are taken into account. Perhaps this issue can be addressed by introducing non-shift symmetric terms in the K or G 3 functions, relying on lower derivative screening mechanisms.
Finally, the question of the stability could be examined in order to carry further the works presented here. Notably, the scalar field perturbations will propagate in an effective metric that differs from the spacetime metric. This metric can be found in [29]. One should check whether this metric is always Lorentzian, and if the associated causal cone is compatible with the causal cone of the spacetime metric, in the sense of refs. [36,37]. Another possible extension of our work would be to consider the thermodynamics of the black hole solutions presented here. The temperature and entropy could be computed along the lines of, e.g., [38]. There are no conceptual difficulties in applying such a procedure; however, since no analytic solution is available in our case, this analysis would require a full numerical sampling over the parameter space of the model, eq. (19), which is beyond the scope of this paper. by an STFC Consolidated Grant.