Exploring the black hole spectrum of axionic Horndeski theory

We provide a ``user-friendly'' algorithm to systematically and rapidly obtain exact planar black hole solutions in the Einstein-Maxwell theory deformed by the most general shift- and reflection-symmetric Horndeski sector where the usual Galileon is replaced by a tuple of scalars with profiles linear in the coordinates of the transverse manifold. Under precise assumptions, these axion backgrounds break the translational invariance of the system, causing momentum dissipation in the holographically dual field theory. The success of the method relies on the simple realization that the bulk equations of motion become more tractable when written in terms of the axions kinetic terms, instead of the radial coordinate. Showcasing this particularly efficient recipe, we derive novel asymptotically AdS black holes, and show that their extremal counterparts always flow to an $\mathrm{AdS}_2\times \mathbb{R}^2$ infrared fixed point. Additionally, we report an interesting family of new asymptotically Lifshitz black hole solutions with $z>1$. Finally, we discuss the DC transport properties of the dual relativistic field theories in view of possible applications to condensed matter systems.


Introduction
With the advent of the Anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [1] and the subsequent advances in "bottom-up" holography [2,3], hairy black holes have become interesting theoretical laboratories offering a fresh insight into the physics of strongly correlated systems, particularly in the context of an AdS/Condensed Matter Theory (AdS/CMT) duality [4]. Perhaps one of the first and most prominent examples might as well be the holographic superconductor (or more precisely, superfluid) [5] of Hartnoll et al.; a black hole in the bulk, which develops charged hair only at low temperatures, is utilized to successfully describe the formation of a charged condensate in the boundary field theory, and the subsequent superfluid phase transition that takes place.
Especially relevant when dealing with such gauge/gravity models, planar black holes, that is, black holes with a Riemann-flat transverse manifold, have their own place in the "black" landscape of gravitational theories. However, precisely because the transverse manifold has vanishing curvature, not all planar solutions describe black holes, and the task of finding those that ultimately do so, is not always a trivial one. In fact, it is a simple exercise in pure General Relativity (GR), to show that any planar solution inevitably describes a naked central singularity, unless a negative cosmological constant is introduced.
In [6], a detailed analysis has been performed in order to determine what kind of matter fields Einstein's theory might support, for a quite general metric ansatz permitting a weak version of Birkhoff's theorem. A plethora of novel solutions is featured there, describing static hairy black holes, dressed with form fields, in the presence of suitably chosen matter actions. A particularly interesting family of solutions contained in [6] is a flat-horizon one in which the matter content is made of multiple k-form fields, the number of which depends on their rank k and the dimension d of the transverse manifold, such that the overall stress tensor is eventually compatible with the isometries of the latter. Using these k-form fields (of rank not necessarily equal to d), the authors proceed with the construction of four-and higher-dimensional asymptotically AdS (AAdS) black holes possessing primary hair.
Most pertinent to this work, are the four-dimensional black holes in [6], dubbed "axionic black holes", dressed with two constant 3-forms H (i) = p dt ∧ dr ∧ dx i , where i = 1, 2. By means of a Hodge duality these solutions can also be constructed using two constant 1forms F (i) = p dx i . These are the field strengths of scalar fields φ (i) = p x i , defined up to a constant, so that any Lagrange 4-form H (i) ∧ * H (j) δ ij is replaced by a sum of kinetic terms, dφ (i) ∧ * dφ (j) δ ij , which is exactly the main building block in this work. Note here, that these configurations exhibit a peculiar feature; altough the horizon is flat, the inclusion of axions leads to a lapse function with a form reminiscent of hyperbolic black holes, due to the emergence of a negative effective curvature scale, proportional to the axionic charge p. In the following sections, we will actually see that such a phenomenon is basically tied to the presence of standard kinetic terms with the "correct" sign in the action.
Solutions with axionic charges become particularly relevant and interesting when taken as gravitational duals for certain strongly coupled quantum field theories where momentum is dissipated. In particular, in recent years, a lot of effort has been devoted to the usage of the gauge-gravity duality to understand condensed matter systems at strong coupling in which quasiparticles are absent and strong correlations are active [4,7,8]. Part of the interest has been motivated by the unusual (e.g. non-Fermi liquids like) transport properties of the so-called strange metals and the origin of superconductivity in their high-T c counterparts. A crucial step in making this tool operative and therefore applicable to realistic condensed matter systems was to realize the fundamental importance of translational symmetry breaking and the consequent dissipation of charge carriers momentum. Abandoning the full diffeomorphism group invariance in the gravitational bulk was the key step in this direction. The explicit breaking of spatial translations of the dual field theory directly implies the presence of a finite graviton mass in the bulk [9], and it elevates massive gravity as the universal low-energy description for condensed matter systems with finite electric conductivity [10]. Interestingly, this statement has been proven exactly, at least at the perturbative level, starting from a proper periodically modulated holographic lattice [11].
A very convenient way to formalize theories of massive gravity is via the so-called Stückelberg mechanism [12,13]. In its simplest incarnation, this formalism makes the new extra degrees of freedom explicit; they directly appear in the action in the form of a set of massless scalars whose expectation values are linearly proportional to the spacetime coordinates. In order to keep the total energy of the boundary field theory conserved, and therefore the time component of the boundary Ward identity unmodified, i.e., ∂ a T at = 0, we do not break time diffeomorphisms in the bulk. This can be ensured by considering a minimal set of scalar fields given by: where the index j = 1, ..., d only spans the coordinates of the transverse manifold, and δ I i is an isomorphism mapping to an internal space of dimension d equipped with a Euclidean metric δ IJ . Notice that these are precisely the fields we previously discussed. We remark here that this choice is not the most general, but the only one retaining the rotational invariance in the spatial sub-manifold 1 . This approach was first introduced in the applied holography community in [16] and later generalized and understood in detail in [10,17]. Nowadays, it is customary to call this setup "holographic axion model" [18], and to recognize it as the most versatile, simple and efficient holographic framework for breaking translational invariance. Importantly, this model is not only able to introduce momentum dissipation, but also the spontaneous (phonons dynamics) [19] and pseudo-spontaneous [20,21] breaking mechanisms.
The big motivation to look into these models is given by the challenge to reproduce and understand the exotic transport properties of strange metals, strongly correlated materials which do not follow standard Fermi Liquids theory. Despite some initial excitement [22,23], the problem still remains unsolved [24]. Additionally, almost all condensed matter systems do not display relativistic invariance and most of them exhibit critical points with nonrelativistic scalings (see Section 2.1.1 in [3]), the so-called Lifshitz scalings: with z > 1 [25]. Having a non-trivial Lifshitz scaling exponent seems to play an important role in the phenomenology of cuprates [26]. Its holographic incarnation is manifested in the so-called Lifshitz spacetime [27], which has opened a completely new avenue for the study of strongly coupled non-relativistic field theories using the holographic duality. The whole dictionary for gravity duals with Lifshitz asymptotics has been recently formalized in [28] and it presents substantial differences with the standard procedure in AAdS bulk spacetimes. Given these facts, it is worth to consider more general and complicated gravitational models. In this direction, Horndeski theory [29] is a promising framework since it constitutes a natural extension of the aforementioned holographic axion model, and it has indeed already been considered before [30][31][32][33][34][35][36][37][38][39].
Horndeski theory is the most general four-dimensional scalar-tensor theory with equations of motion up to second order in derivatives of the fields. Although the action contains higher-derivative self-interactions of the scalar which in general contribute higher-order terms to the field equations, it also contains gravitational counterterms, the contribution of which results in finally having a healthy set of second-order differential equations. The shift-symmetric Horndeski action can be written as: in its modern generalized Galileon formulation [40][41][42], with Here, R is the Ricci scalar, G µν is the Einstein tensor, and G nX ≡ ∂G n /∂X is a convenient notation. The invariance of (1.3) under constant shifts of the scalar field, φ → φ + c, requires that the usual dependence of the Horndeski functions G n on φ drops out; these functions only depend on the kinetic term X ≡ −(∂φ) 2 /2, and thus, φ is massless. In this work, we wish to explore static planar black holes with scalar backgrounds of the form (1.1) in the framework of Horndeski gravity. To do so, we have to rewrite (1.3) so that it may support such solutions in the first place, i.e., we must promote the Horndeski scalar φ to a doublet of scalar fields {φ I | I = 1, 2}. Doing so, we further wish to secure the invariance of the action under both spacetime and internal transformations which means that the action itself should be a true coordinate and internal scalar with no free indices. Internal transformations here will actually be appropriate scalar field reparametrizations, as we will later see in detail. Clearly, this can only happen if we discard L 3 and L 5 , or equivalently, in terms of the original Horndeski action (1.3), if we only keep the sector invariant under reflection transformations φ → −φ. The explicit form of the truncated action for the doublet will be given in the next section, and we will refer to it as axionic Horndeski theory. Note here that one could judiciously question the necessity of forming internal scalars in order to obtain bulk solutions with scalar backgrounds as in (1.1). Indeed, we might as well introduce two scalar fields φ I by just adding two copies of (1.3), one for each scalar 2 . However, it is interesting, that even in this case, in order for (1.3) to support solutions with scalar backgrounds linear in the coordinates of the transverse manifold, one needs to perform a significant truncation, namely, L 3 and L 5 must again be discarded, but for a different reason now. The reason has to do with the equations of motion for the scalar fields, which are just on-shell closure statements for the two Noether currents associated with the aforementioned shift symmetry. By direct inspection, one can verify that for a profile (1.1), the equation of motion for the Ith scalar field cannot be satisfied essentially due to the presence of a radially dependent radial component of the Ith current. It can be further shown that this radial component is solely a variational artefact of having L 3 and L 5 in the action; once we discard those, the Ith current does only possess a radially dependent component in the Ith direction of the transverse manifold, and therefore, the corresponding scalar equation of motion is identically satisfied. Nevertheless, we won't be dealing with such a prescription, and we just mention the above for completeness.
Finally, the main goal of this work is to establish a systematic approach to hairy planar black holes in axionic Horndeski theory. In simpler words, if such exact solutions exist, we wish to be able to derive them with minimal effort by just specifying the Horndeski functions, and executing some straightforward process steps. In this sense, we want to construct a step sequence which results in solutions once executed, a solution-generating algorithm if you wish. Note that although we only display certain examples, this algorithm can in general produce many more exact solutions with the aforementioned characteristics, and in that fashion, the principal result of this manuscript is a rather general one. Based on the logic behind this prescription, and armed only with a minimal set of initial assumptions, we will also be able to pull off a complete derivation of the DC transport matrix for relativistic field theory (holographic) duals, without the need to fix the G n 's a priori. Although this part will not be fully developed in terms of an exhaustive analysis, this is nevertheless a potentially useful result, as universal linear-response features might possibly be extracted from it.
Plan of this work.
The manuscript is organized as follows: in section 2 we present our general setup and its main features; in section 3 we present the novel algorithm which can be used to generate planar black hole solutions with axionic hair, this being the main result of our work; in section 4 we provide two examples of how this algorithm works; in section 5 we initiate the study of the relativistic field theory duals of these solutions by deriving their DC transport properties; finally, in section 6 we conclude and discuss possible future directions.

The general setup: theory and field equations
In this work, we deform the Einstein-Maxwell (EM) action by adding to it a Horndeski action with reflection and shift symmetry. Instead of a single scalar field, we introduce dmany axionic scalar fields, or simply axions, in general d + 2 spacetime dimensions, packed in a d-tuple φ = (φ 1 , ..., φ d ) which transforms as a vector under the action of an internal group O(d). This "flavor" vector can be used to construct a scalar quantityX: where I IJ ≡ − 1 2 g µν ∂ µ φ I ∂ ν φ J , capital Latin indices I, J, ... are internal-space indices ranging from 1 to d, and the trace is taken with respect to the reference metric δ IJ = (1, 1, ..., 1) IJ , which is the canonical static metric in d-dimensional Euclidean internal space. 3 This is basically the main building block of our theory, as all Horndeski functions will be functions of (2.1). Notice that for d > 1, this is not the most exhaustive choice; already for d = 2, there is another algebraically independent contraction: det I = 1 2 Tr I Tr I − Tr I 2 . (2.2) 3 Einstein's summation convention is adopted also for internal-space indices.
Indeed, for a d-dimensional transverse manifold, I is a d × d matrix, and Tr I n with n ≤ d are independent scalar objects which can, and in general have to, appear in the action via appropriate combinations. 4 Nevertheless, in this work we do not consider other combinations.
More concretely, we choose to work with the following action principle Here, we have set the coupling in front of the Einstein-Hilbert piece in (2.4) to unity which amounts to setting M d P l = 16π, together with natural units. Moreover, Λ denotes the bulk cosmological constant, and the tensor F µν = ∇ µ A ν − ∇ ν A µ stands for the field strength of a bulk U (1) gauge field A µ . We also use the convenient notation: where the internal index should always precede the coordinate ones. Finally, the last piece in (2.5) can be cast into the more familiar form: and one may then immediately recognize the quartic sector of (1.3). Let us first discuss the symmetries of (2.3). Apart from being a coordinate scalar, the action is also an internal-space scalar. Therefore, it is invariant under general coordinate transformations, but also under scalar-field reparametrizations which are isometries of the Euclidean plane. The latter are transformations where Λ is a static matrix in O(d) and c is a constant shift vector. As mentioned in the introduction, the action (2.5) is the only sector of the (d + 2)-dimensional version of (1.3), which, after promoting the Galileon to a d-tuple, enjoys these symmetries. Indeed, in the other sectors, φ appears an odd number of times which implies that these coordinate scalars would transform as O(d) vectors once we replace the scalar field by the flavor vector. For example, the higher-derivative self-interaction terms in L 5 would in general have three free internal indices, but only two of them could be contracted with the reference metric, thus leaving us with one free internal index. Additionally, note here that once a background configuration is chosen for the scalar fields, the full symmetry will be broken down to a diagonal subgroup of internal and spacetime symmetries (see [10] for more details).
To extract the bulk equations of motion, we vary (2.3) with respect to the metric, the gauge field and the axions φ I . Variation with respect to the axions, using represents the Ith Noether current associated with the constant shift φ I → φ I + c I . Hence, the equations of motion for the d-many axions are just (on-shell) closure statements for the d-many such currents. Moving on, the gauge-field equations of motion are the standard ones: Finally, the modified Einstein equations, obtained by varying (2.3) with respect to the metric, are To get the above, we also used the fact that for any scalar function G(X), it holds that (2.14)

A solution-generating algorithm
The main task of this work is to construct a systematic way of obtaining all solutions following from (2.3), which further satisfy the following requirements: (i) The bulk geometry is four-dimensional, represented by the static ansatz: where h, f are to be found by solving the bulk equations of motion. The coordinates x i of the transverse manifold are considered noncompact, and in this regard, the spacetime is truly plane-symmetric. The boundary is located at radial infinity, r = ∞.
(ii) The gauge potential is a function of only the radial coordinate r with the only nonvanishing component being the temporal one, A t (r) = 0. In particular, the gauge field reads A µ = (a(r), 0, 0, 0).
(iii) As already anticipated in (1.1), the bulk solution for the axions is With these requirements, one can also straightforwardly verify that E µν is compatible with the symmetries of (3.1), i.e., £ ξ E µν = 0 where the Lie derivative is taken with respect to all the corresponding Killing vectors. Note, however, that because of (1.1), £ ξ φ I = 0, and in this sense, the full solution (that is, including the gauge and scalar backgrounds) is only static. Using the algorithm, we will be able to generate plane-symmetric solutions to the bulk equations of motion, but we won't be able to tell if these describe planar black holes. This is an issue that has to be treated separately for each case. We will put the algorithm to use, providing novel closed-form hairy black holes which we classify according to the form of G 4 . We will refer to these black hole solutions as axionic black branes [46]. Moreover, in this work we will not make statements about the physical permissibility of these configurations in terms of instabilities, energy conditions etc. In this vast landscape of solutions, it is reasonable to expect the existence of pathological setups which should be avoided. Given these premises, we havē where ′ ≡ ∂/∂r . Here, we can make a first observation which will prove to be crucial for the success of the algorithm. The first equation in (3.2) can be inverted in order to find withX < 0 by definition. As a consequence of the reflection symmetry, p will always appear inside a modulus, and henceforth, we take p > 0 without loss of generality. It turns out that if we treatX as a variable, instead of the radial coordinate r, the equations of motion become much more tractable, and an algorithmic approach can be established. Radial derivatives of any function G(r) can be written asX derivatives, e.g., and we will use this fact, together with (3.3), to transform ordinary differential equations (ODEs) for r to ODEs forX. Let us first start with the equations of motion for the axions. In the previous section, we saw that the Ith axion must satisfy ∂ µ ( √ −gJ µ I ) = 0. To verify that φ I = p δ I i x i is a solution, we repeat here the argument made in the introduction of this work. Evaluating (2.10), one finds that where j(r) is some complicated function of the radial coordinate. Therefore, it immediately follows that identically. Next, the equations of motion for the bulk gauge field can be solved for a closed-form expression of the electric field: where the integration constant q e is related to the charge density in the dual field theory.
We note here that one may also include a magnetic charge by considering a dyonic profile for A µ . Since the electromagnetic sector of the action consists of just the standard Maxwell term, at least at the level of the background solutions, the presence of magnetic components results only in an additive contribution; that is to say, solutions with a dyon can be obtained from the electrically charged solutions we present, by sending q 2 e → q 2 e + q 2 m . As expected, dyonic solutions will be invariant under the duality exchange q e ↔ q m . Although this may be not that interesting from a gravity point of view, considering dyons opens several new possibilities for the transport properties of the dual field theory. We leave the study of the dyonic solutions for the near future.
Proceeding with the modified Einstein equations (2.12), it is possible to show that, given the solutions φ I = p δ I i x i and (3.7), the ii components are not independent, and therefore we only need to satisfy E tt = 0 = E rr in order to solve the full system. This has been proven in full generality in [47][48][49], where the authors show how to generalize a Schwarzschild-like ansatz in order to construct static black hole solutions for any Lagrangian density (describing a metric theory of gravity) being a functional of the Riemann tensor and its derivatives of any order, arbitrary k-form fields and scalar fields. The metric ansatz is generalized by considering an isotropy-irreducible homogeneous base space which allows, before specifying the particular theory under scrutiny, to ensure that solutions of the system can be found by solving two ODEs, one for g tt , and the other for g rr (please refer to [49] for further details).
The tt component yields a first-order inhomogeneous ODE for f with function coefficients α n . In theX formalism, it reads Then, E rr yields another ODE, for h this time, which acquires the particularly simple form The general solution to (3.8) is obtained via the integrating-factor method: provided that α 2 = 0, whereas (3.12) can be solved for Here, both f 0 and h 0 are integration constants.
To simplify things, we now set h(r) = C(r)f (r), and observe that a master equation for C can be found. It is given by Remarkably, C knows nothing of k-essence, and it can be directly integrated from the above equation once G 4 is chosen. Of course, a closed form result depends on whether theX integral of the right-hand side (rhs) of (3.16) has a closed form expression, or not. Note here that the denominator of (3.16) is ∝ α 2 , and thus, nonzero. Another interesting remark can be made. To obtain solutions with g rr = −g tt , the rhs of the above equation must vanish, and this can only happen if G 4 is constant, or for G 4 ∝ √ −X. The general black brane solution for this particular model has been already found in previous work [31], but with a slightly different action and building blocks. 5 Now that we know C, we can find a better form for f by massaging (3.14) in order to arrive at The above form is much more tractable now, and we will see that with (3.17) and (3.16) at hand, we can rapidly obtain solutions which would otherwise be hard to directly integrate from the field equations.
Let us now assume that the solution describes a black brane with its event horizon located at r = r 0 . Via the standard Wick-rotation method and the avoidance of the conical singularity in the Euclidean section, we can find the Hawking temperature which in theX formalism assumes the form: withX 0 ≡X(r 0 ). In a less obscure format, in terms of r 0 , where all functions are evaluated at r 0 . Notice that C(r 0 ) must be positive for T H to be real; needless to say, we must choose C(r) > 0 for all r ≥ r 0 in order to have positive −g tt , g rr in the exterior. At this early stage, we can also extract an expression for the Wald entropy [50]: whereǫ µν is the unit binormal to the bifurcation surface Σ with normalizationǫ µνǫ µν = −2, and we see that the result is completely ignorant of k-essence since the latter couples minimally to gravity. Note that one can define an effective running gravitational constant, in particular, thereby re-expressing (3.23) as In units G ef f (r 0 ) = 1 the renown "quarter law" is then manifest [51]. However, this definition must come with G 4 > 0, otherwise G ef f will become negative at some r. Although we will not define such an effective object, we will take G 4 to be positive in order to avoid imposing more delicate conditions on r 0 . Finally, bear in mind that the plane extends to infinity, and thus, for −ℓ < x i < ℓ with ℓ → ∞, The quantity of interest will therefore be the entropy density: and by "finite entropy" in the following text we will mean a finites.
Since we know the Hawking temperature, we can take things one step further, and assume the existence of extremal black holes, zero-temperature black hole solutions to the bulk equations of motion. Given a bulk solution at Hawking temperature T H , its extremal version, if it exists, can be found by solving T H = 0 for r e , the largest positive real root of T H which also is a double zero of h and f . In equivalent words, we need to find the largest solution to − q 2 e − 4Λr 4 + 2Kr 4 = 0. (3.28) It is interesting that G 4 does not enter in the above, and we can conclude that extremal branes come with sizes independent of the quartic-sector contributions, although the latter are expected to affect the mass of the extremal configurations. Before closing this section, we would like to see what information we can get from the zero-temperature Einstein equations without knowing the explicit forms of h, f that solve the field equations in the whole bulk.
To that end, let us work with the near-horizon approximation f (r) = f ′′ (r e )(ǫz) 2 /2 + O(ǫ 3 ) and similarly for h, where we defined ǫz ≡ r − r e . Given an electric field (3.7) and the axion backgrounds (1.1), we do a series expansion of the metric-field equations about ǫ = 0 with z held fixed. To leading order we find that the tt and rr components are simply the statement T H (r e ) = 0. However, nontrivial constraints are coming from the ii components, in particular, 29) all functions evaluated at r e . Additionally, (3.1) becomes where h ′′ (r e ) = C(r e )f ′′ (r e ). We now need to define a new time coordinate τ in order to absorb ǫ so that we can then send ǫ → 0 with τ, z held fixed obtaining the extremal IR metric. Since C(r e ) is ultimately a number, it can also be absorbed into τ by t → τ /( C(r e )ǫ). Combining this with the optional reparametrization u(z) ≡ 2/(f ′′ (r e )z), leads us to the near-horizon AdS 2 × R 2 form: where we have identified . (3.32) The above makes sense for f ′′ (r e ) > 0. Note that at this stage, this is only a solution in the deep infra-red. However, once we find a bulk solution with a zero-temperature configuration, the extremal setup will automatically describe a geometry interpolating between (3.31) in the far IR and some solution-dependent asymptotic geometry in the UV. 6 Finally notice from (3.29) that the extremal geometry exists also in the absence of an electric charge, purely supported by the axionic parameter p, provided that KX = 0.
The "user-friendly" algorithm To conclude, let us give a "user-friendly" version of what we just presented with the visual aid of the flowchart in figure 1. Given the solutions (3.7) and φ I = p δ I i x i , we can group the total process of solving the modified Einstein equations into three phases: (ii) Integrate the rhs of (3.16) with respect toX, then exponentiate to obtain C.
This results in an expression for C.
(ii) Evaluate the integral in the rhs of (3.18) to obtain k.
This results in an expression for k.  Note here that phase 2 can be skipped, and the solution can be given in terms of k. In what follows, we will initially skip phase 2 in order to do a preliminary analysis of the solution, but we will always complete this phase afterwards. In this sense, solutions should be first classified with respect to G 4 .
Further data can be also collected by direct evaluation: • Choose G 4 and K, complete phase 1, and evaluate the rhs of (3.20) to obtain T H . 6 UV and IR do of course make sense only in the language of the dual field theory. That being the case, the domain wall picture is, roughly speaking, a "geometrization" of the RG flow in the dual system.  Figure 1: Flowchart description of the algorithm. We can group the steps into three phases which are indicated by numbers enclosed in disks. Ellipses indicate start/stop, trapeziums are for manual input, rectangles denote process steps, parallelograms are for output and the circle is a combination of outputs from phases 1 and 2.
• Choose G 4 and evaluate the rhs of (3.23) to obtain s W .
• Choose K and solve (3.28) for r e (if applicable).
• Choose K and G 4 , determine r e from previous process, evaluate the rhs of (3.29), and substitute it into (3.32) to find the size of the fixed-point AdS 2 (if applicable).
In the next section, we initiate a first exploration of the vast landscape of models captured by (2.3), with the assistance of the algorithm.

The G 4 ∝X class
As our first example, we set G 4 = γX. Completing phases 1 and 3, we quickly arrive at the full solution: Here, µ is the chemical potential in the dual field theory, given in terms of horizon data so that a(r 0 ) = 0, i.e., µ = π 2γ q e erfi γ 2 The integration constant M is related to f 0 in (3.17) via M = pf 0 , and it is also related to the physical regulated mass per unit area in the (x 1 , x 2 ) plane. Additionally, we have also definedγ ≡ γp 2 , and D denotes the Dawson function [52]: We take G 4 > 0 which amounts to γ < 0. This is consistent with the reality of (4.2), and it further guarantees a non-negative Wald entropy.
Expanding at large r, we have where k r→∞ denotes the full asymptotic series expansion of k with leading-order term k r→∞ . Note that r → ∞ amounts toX → 0. Let us assume that k admits a power series expansion about the boundary with k (0) r→∞ ∝ r ∆ . To leading power, the first terms in (4.6) and (4.7) grow as ∼ r ∆−1 . If ∆ < 3 the asymptotic geometry is locally AdS once we fix Λ = −3, with metric: 7 in terms of the holographic radial coordinate u(r) ≡ r −1 . If ∆ = 3, then we observe that and one can still achieve the asymptotic form (4.8) provided Λ ef f < 0. On the other hand, if ∆ > 3, we have that g rr ∼ r ∆−1 close to the boundary. In [53], it was shown that this asymptotic behavior is tied to diverging curvature invariants as r → ∞. This can be seen by writing the Riemann tensor in the static orthonormal frame, the diagonal tetrad, and asymptotically expanding it, only to find out that positive powers of the radial coordinate appear in its components. Consequently, any invariant built from the Riemann tensor will unavoidably diverge as r → ∞, rendering the solution pathological. Especially from a holographic perspective, this is ultimately unwanted since gravitation would not become decoupled at the boundary. Therefore, one should only consider solutions with ∆ ≤ 3.
So far, Eq. (4.2) is a solution in models which may or may not contain kinetic terms for the axions, collectively denoted byX. To find the contribution of kinetic terms to the solution, we set K =X +K, and perform phase 2 in order to obtain: and thus, the presence of kinetic terms in (2.5) corresponds to a contribution at the level of (4.2), which asymptotically is in O(1). We mention here that the highly relevant and well-studied solution in [36] follows from (4.2) by choosing K ∝X, or equivalently,K = 0. Hence, we have provided a straightforward extension of the solutions in [36] for literally any reasonable functionK with contributionk given by (4.11).
There are many choices that yield a closed-form expression for the integral in the rhs of (4.11), and as an example, we will chooseK = β(−X) n which leads tõ Here, Γ is the incomplete Gamma function [52] with the series expansion (4.14) Remember that we previously demanded ∆ < 3. This amounts to n > 0, and we have three asymptotic branches fork, namely, where the proportionality factor is positive. Since n > 0, (4.2) can only describe a family of AAdS black branes when Λ < 0, the horizon structure of which we are going to investigate in the remainder of this subsection. Let us set Λ = −3 and obtain the Hawking temperature by evaluating (3.20): The first step will be to determine the location of the extremal horizon by finding the largest solution to 12r 4 + 2βp 2n r 2(2−n) − 2p 2 r 2 − q 2 e = 0 (4.17) This function will come in handy in the analysis that follows. We study (4.17) numerically, finding cases with even two positive real roots (see figure 2). In general for l distinct positive real roots, the largest of them will be the location of the event horizon of the extremal black brane, whereas the remaining (l − 1)-many ones are expected to give rise to states with up to l-many inner horizons. Let us explicitly show this for the case at hand by using (4.18). Assume there are two roots, r 1 and r 2 , with r 2 > r 1 . Both are roots of T H , and thus, are also double roots of h, f . Clearly, r e = r 2 . Additionally, they correspond to extrema of M. The latter goes to ∞ as r → ∞ by definition, while it goes to a finite valueM as r → 0, in particular,M The extremal black brane has M e ≡ M(r e ), with a naked singularity for all M < M e . However, f goes as Since M has two extrema for r > 0, and M e is a global minimum, it has to be that M 1 is a local maximum, and thus, M e <M < M 1 . This information is enough to conclude that: (i) For M < M e , the solution describes a naked singularity.
(ii) For M = M e , the solution describes an extremal black brane with its event horizon located at r = r e .
(iii) For M e < M ≤M, the solution describes a black brane which further develops an inner horizon. The causal structure is similar to that of Reissner-Nordström black holes.
(iv) ForM < M < M 1 , the solution describes a black brane with two inner horizons. This is the picture of a "black hole inside a black hole" [54].
(v) For M = M 1 , the two inner horizons coalesce into one at r = r 1 .
(vi) For M > M 1 , the solution describes a black brane with a regular event horizon.
In the present work, we have not studied in detail the nature of the geometry inside the black hole horizon, but it would be interesting to perform such analysis on the lines of [55][56][57][58]. It is reasonable to expect that the Horndeski terms play a role in such a game as already emphasized in [59]. All of these extremal black branes are finite-entropy domain wall solutions interpolating between a unit-radius AdS 4 in the UV and an AdS 2 × R 2 IR geometry, where the size of AdS 2 can be numerically obtained for all cases via (3.32). An exact example of this may be useful. Let us set n = 2. For Then, using (3.32), we find that The fact that there is no dependence on γ is consistent with (3.29), where the denominator simply becomes r 4 e . As a concluding remark, we mention here that we have performed a surfacial analysis of a very specific family in the G 4 ∝X class, given by K =X + β(−X) n . There are many other choices K =X +K which still yield exact solutions, and it would be interesting to investigate them.

The G 4 ∝ (−X) m class
As our second example, we set G 4 = γ(−X) m . Completing phases 1 and 3 in figure 1, we immediately arrive at the full bulk solution: Here c 0 is an integration constant related to the chemical potential µ in the dual field theory, while The integration constant M is related to f 0 in (3.17) via M = pf 0 , and again, it is also related to the mass (per unit area) of the solution. The symbol F 2,1 denotes Euler's hypergeometric function [52]: with On the other hand, if z → ∞, one needs to perform Pfaff's transformation [52] in order to send z → z −1 , provided that the phase of −z is less than π. Then, The solution is well-defined for C > 0 (which amounts to c m < 0) and all m such that the hypergeometric functions are well-defined. In particular, the former restriction takes the form m < 1 under the entropy requirement G 4 ≥ 0 (meaning γ>0), whereas the latter restriction can be shown to be Here, we used ∧, ∨ as standard symbols for logical operations "and","or", respectively. This condition looks absolutely horrible and cryptic, but bear with us since will eventually visualize it in just a bit. Now, z = c m in the above expansions, and c m → 0 as r → ∞ for m > 0. Therefore, to leading order at large r, and we may use the arguments of the previous subsection regarding the Riemann tensor to constrain the leading mode in k r→∞ . Practically, this means that for any positive 0 < m < 1, we only consider AAdS solutions. Now, for m < 0, c m → ∞ as r → ∞, and the picture is quite different: thus, the only asymptotic geometry of interest (with constant curvature invariants at radial infinity) will be Lifshitz. However, before we start fixing parameters, let us completely determine the solution by choosing K =X + β(−X) n and completing the missing phase 2 in figure 1. Doing so, we find the contribution: with c m , v, w defined in (4.25). This particular form implies the identification: For (4.32) to be well-defined, we need to impose further restrictions which we do not write down explicitly; basically, we take the union of these conditions with (4.29), and we depict the logical negation of that union as white space in figure 3, together with exclusions coming from filtering out the cases with diverging tidal forces close to the boundary. Most importantly, in figure 3 we display the three distinct kinds of asymptotic forms we end up with. Notice that we have severly restricted the (m, n) parameter space. As previously anticipated, for negative m there exists a family of asymptotic Lifshitz (ALif) geometries. We will show that this corresponds to the lower boundary of the blue region in figure 3. Lifshitz spacetimes are described by the metric [60]: where z is the dynamical critical exponent, and we have set L = 1. The key-feature of (4.34) is its symmetry under the anisotropic scale transformations Here, we expect our black branes to be locally asymptotically of the form (4.34) with z being a function of m, n, under certain conditions. First, if h is going to grow faster than f , it has to be m < 0. Then, if f is to grow as ∼ r 2 , it further has to be that n < 0 and m = n. As we have set L = 1, the coefficient of the g rr leading mode has to be set to unity which implies β = 2γ(3 − 6m + 4m 2 ). (4.36) This ensures that g rr ∼ r 2 (1 + ...) at large r, where the dots denote terms in orders of growth less than O(1). Having secured the desired growth for f , we can identify with z > 1 amounting to m < 0, which is exactly the branch we are studying. Additionally, using (4.36) and setting we can make sure that −g tt ∼ r 2z (1 + ...) as r → ∞. Summarizing, the tuning: together with the identification: which is (4.37) inverted, give rise to axionic black branes with a metric locally acquiring the asymptotic form (4.34). To the best of our knowledge, this is the only exact Lifshitz-like family in Horndeski gravity with arbitrary z > 1, 8 where z is fixed once the degree of the G 4 monomial is chosen. Now, we turn our attention to (deep) IR geometries which can act as ground states. One can show that for (4.24), it is impossible to construct a regular core near the origin of the radial coordinate; the curvature invariants exhibit poles which cannot be removed by tuning. Hence, the only candidate left is the extremal IR geometry (3.31). As we discussed in the previous section, the location at which the extremal black branes develop a horizon, is the same for all G 4 classes with a common choice for K, a fact which saves us a lot of time, since we have already performed this analysis in the last subsection where we also chose K =X + β(−X) n . In the current case, it makes sense to extend figure 2 to the negative n half plane which is occupied by configurations where (4.17) has one positive real solution corresponding to the location of the extremal horizon. Apart from the usual AdS 4 → AdS 2 × R 2 flow for 1 > m > 0, n > 0, here we also have extremal finite-entropy bulk geometries interpolating between a Lifshitz asymptotic geometry and AdS 2 × R 2 in the infra-red. For these zero-temperature configurations, we numerically determine the extremal value of M , M e (z), plotting it in figure 4. Close to z = 1, the z discontinuities become more frequent, and M e exhibits a peculiar oscillating behavior. As z grows, M e goes to some constant value. These wild fluctuations also cease to exist once we go to the electrically neutral version of these extremal branes, as can be seen in figure 4. Overall, the charge seems to have an uplifting role, whereas a positive Λ seems to have the opposite effect, pushing M e to lower values. Note that Λ ≥ 0 is allowed in this case, since the bulk cosmological constant is unrelated to the locally asymptotically Lifshitz form of the metric,  Figure 5: Plotting −g tt (orange) and g rr (blue) versus the radial coordinate r for neutral ALif extremal black branes with Λ = 0 at p = 1. From high to low opacity, z = 3/2, z = 2 and z = 5/2, respectively. a fact made obvious by its absence in the tuning (4.39). Also, negative values for M e , and thus, negative values for M > M e , is a known feature for planar solutions with axionic charges (see for example [6]).
Let us close this subsection with the easiest exact example, the case with q e = 0 = Λ. The extremal ALif black brane will develop a horizon at (4.41) the latter supported by the axion flux p. This is the exact solution to (4.17). The nearhorizon metric will have the form (3.31), with the size of the IR AdS 2 given by obtained by evaluating (3.32), 9 where we expressed the result in terms of m(z) (see (4.40)). The behavior of the two metric functions, h and f , can be seen in figure 5. The full (extremal) solution is then interpolating between a UV geometry (4.34) with the Lifshitz scaling (4.35), at least up to some appropriate cut-off, and an IR geometry (3.31) with the radius (of curvature) of AdS 2 given by the square root of (4.42), while its entropy is also finite:s As a concluding statement, we also remark here that (4.24) matches the solution in [31] for m = 1/2 and K =X, at least up to some insignificant arithmetic differences, which are due to the slightly different action principle used there.

Holographic DC conductivities
In this section, we consider a general black-brane background like those in the previous section with the further requirement that the UV fixed point of the dual field theory is relativistic invariant, i.e., the bulk geometry is (locally) asymptotically anti-de Sitter. We assume the presence of an event horizon at some radius r = r 0 . We use an Eddington-Finkelstein (EF) coordinate system in which the regularity of the metric (3.1) in the vicinity of the horizon is manifest. The Hawking temperature is given by

2)
A finite chemical potential, defined as is introduced in the dual field theory. In the rest of this section, we will nevertheless work in the thermodynamic ensemble where the charge density is kept fixed, instead of the chemical potential. The electric charge density is given by where in the last equation we used (3.7).
We also once again write h = Cf . Looking at the Maxwell equations of motion, at linear level there is only one nontrivial component, the one in the x 1 direction: where, using the background equations, Here, J x 1 represents the time-independent electric current in the x 1 direction when evaluated at r → ∞. Notice that time dependence will drop out from (5.10), if We will keep this in mind. Turning our attention to the linearized Einstein equations, we have two nontrivial components, tx 1 and rx 1 . Considering E tx 1 first, and by using the (linearized) Maxwell equations of motion (5.10), we can show that where the γ n 's are complicated functions of r. Their explicit form is irrelevant to our purpose for the moment. The key step here is to manifestly write this equation as ∂ rQ = 0 which amounts toQ being r-independent.
The concrete expression forQ turns out to be: where U and V are to be determined by subtracting (5.12) from the radial derivative of (5.13) and solving the result order by order in derivatives of δg tx 1 . Doing so, we first find that the vanishing of the coefficient in front of δg ′′ tx 1 forces the relation (5.14) Using this information, we can deduce that the vanishing of the coefficient in front of δg ′ tx 1 can only be satisfied if Since the γ n 's contain Horndeski functions, it will be convenient to write everything in terms ofX, and only then integrate, finding the solution This also guarantees that the coefficient in front of δg tx 1 vanishes, and results iñ where is the time-independent piece of the thermal current in the x 1 direction when evaluated at the boundary. Also taking into account (5.11), the time dependence ofQ completely drops out if we choose We set c 1 = −E with E parametrizing an electric field deformation. On the other hand, ζ parametrizes a time-dependent heat source, a thermal gradient ∼ ∇T . Henceforth, g 1 , g 2 as in (5.19).
We also need to solve the rx 1 component of the linearized Einstein equations. This is an expression algebraic in the rx 1 mode which has the solution: Observe that at r = r 0 , we have all functions evaluated at r 0 . Clearly, these are finite results, meaning that H rx 1 will diverge as ∼ f −1 when r → r 0 , provided that χ is smooth analytic there. Note also that (5.20) solves the equation of motion for χ.
We must now ensure that all perturbations are well-behaved in the bulk, which means that certain boundary conditions have to be imposed, both at the horizon of the black brane, as well as at infinity. Starting with boundary conditions as r → ∞, from the gauge equations of motion, we see that a sufficient fall-off for a x 1 is ∼ J x 1 r −1 provided that H tx 1 does not yield a non-normalizable mode. The last requirement ensures that only the timedependent piece of δg tx 1 acts as a holographic source for the heat current at r → ∞. Close to the conformal boundary, E tx 1 = 0 has two independent solutions, yielding where c ± are constants of integration, w ± are some functions of γ, ∆, and we assumed that G 4 can be expanded in power series about infinity with leading mode ∝ γr ∆ , together with f, h ∼ r 2 , in order to get the above result. For ∆ < 0, we can expand the hypergeometric function using (4.27); we get a non-normalizable mode ∝ c + r 2 which we kill by setting c + = 0, and a normalizable mode with fall-off in O(r −1 ). This behavior associates with a G 4 which goes to infinity as some positive power ofX. Also, for (5.25) to be well defined in this case, γ(∆ + 2) ≥ 0. Assuming γ is a coupling constant in front of the quartic sector, we saw that the non-negative entropy condition suggests that its sign is fixed positive. 10 This means that −2 < ∆ < 0. But what about ∆ > 0? Assuming K also admits a power series expansion about the boundary, led by ∝ βr Θ , we find that satisfying the bulk equations of motion at an AdS 4 UV fixed point, amounts to ∆, Θ < 0. Therefore, we cannot have ∆ > 0 and AdS asymptopia together, a fact which ensures that the fall-off of H tx 1 is indeed ∼ r −1 . It is satisfactory to verify (as a cross-check) that the asymptotic power restrictions we found here are consistent with those found for the AAdS black branes in section 4. No additional constraints are coming from the asymptotic behavior of δg rx 1 : It can fall off as fast as we wish, with the non-normalizable mode of χ vanishing and its normalizable piece in some arbitrary order of growth strictly less than O(r −1 ). It remains to discuss the behavior of the perturbations at the horizon r = r 0 . To do so we switch to EF coordinates (v, r) with v given by (5.1). Considering the gauge-field perturbation first, we have that Expanding about the horizon, we see that for δA x 1 to be regular there. A valid near-horizon approximation is then Moreover, for δg tx 1 to be regular at the black brane horizon, it must be that where U is some analytic function at r = r 0 with assuming χ is constant there. Note that in order to find this result, we used the approximation (5.29) and the background equations of motion. In this way, the divergence in (5.20) is also remedied, and we finally have a set of well-behaved perturbations in the whole bulk. We can proceed with evaluating J x 1 and Q x 1 at the horizon radius, denoting the values as J and Q , respectively: where all appearing functions are evaluated at r = r 0 . From the above expressions we can read off the components of the holographic DC transport matrix: where ∂ E ≡ ∂/∂E and ∂ ζ ≡ ∂/∂ζ . At temperature T = T H , the electric DC conductivity σ, the thermoelectric conductivities α,ᾱ and the heat conductivity (at zero electric field) κ, are given by 11 respectively. Notice that at zero temperature we have which implies the presence of a resistivity residue. Note that the denominator of the coherent part of σ is related to the (effective) graviton mass evaluated at r = r 0 [9,17,64], thus it better be that Z(r 0 ) > 0, from which it immediately follows that KX > 0 at zero temperature. Moreover, for vanishing heat flows Q = 0, only the incoherent part of σ survives, and thus the conductivity is constant, equal to unity. This is made obvious by the fact that Additionally, the heat conductivity at zero electric current, κ, is given by it can be found by acting with a ζ derivative on Importantly, the discussion just presented is valid only when the scalar operators dual to the axionic fields φ I break the translational invariance of the dual field theory explicitly, causing momentum dissipation. This is in 1-to-1 correspondence with saying that the bulk solution φ I ∼ x I plays the role of an external source for the dual operator. As already explored in the literature [19], this is not necessarily the case. In the opposite scenario, where translations are spontaneously broken, the corresponding DC conductivities would be clearly infinite and must be treated in a separate way not analyzed here.
As an example, let us focus on a specific family of models given by G 4 = γ(−X) m and K =X + β(−X) n . Remember that for these choices, we found AAdS black branes in section 4, provided that 0 < m < 1 and n, γ > 0. We are particularly interested in configurations supporting an extremal profile, meaning that we should avoid the blue region in figure 2. 12 Here, the crucial function Z takes the form: If β < 0, it is guaranteed that Z is positive. For β > 0, one first needs to wisely choose n (based on the value of β) so that the extremal black brane exists. Then, a further condition needs to be imposed, namely, The very large temperature regime is unfortunately not interesting since our model is conformal therein, and it cannot take into account the granularity of matter, for example, the presence of a well-defined lattice spacing. The electric conductivity displays a lower bound, i.e., σ ≥ 1, which holds for all theories with bulk electrodynamics dictated by a pure Maxwell term [65] but that can be easily violated otherwise [66,67]. Notice that at large temperatures (T being the largest scale in the system compared to q e , p, β, γ) the DC conductivity is universally σ = 1. This is just a trivial outcome of the fact that the conformal UV fixed point is reached. It indeed simply coincides with the value one would obtain in a Schwarzschild-AdS background. On the other hand, for small and/or intermediate temperatures, we find a more interesting picture. First of all, the temperature is a strictly increasing function of r 0 , that is to say, it is a positive bijective function for all r 0 > r e which can in theory be inverted to give r 0 (T ). However, for arbitrary m, n we are unable to provide a closed-form expression of the latter, so we will work with parametric plots.
A numerical investigation (see figure 6) reveals that the DC conductivity can be classified into three distinct categories based on the behavior at low/intermediate T /p. These are characterized by: (i) a persistent metallic phase in which dσ/dT < 0 for any T , (ii) an insulator-metal crossover (IMC) where at small temperatures the DC conductivity exhibits a maximum (similarly to what some of us found originally in [17]) and (iii) a metal-insulator crossover (MIC) at low T , followed by an IMC at some intermediate larger scale T /p ∼ 1.
We stress out that one should take the word "insulator" with a grain of salt here. Using the strict condensed matter definition, an insulator is a system for which at T = 0 the T/p σ Figure 6: Electric DC conductivity σ as a function of the dimensionless temperature T /p for q e /p 2 = 1, γ = 0.14, m = 0.13 and: β = 0.78, n = 0.8 (orange -class (iii)), β = −0.2, n = 4 (green -class (ii)), β = −0.5, n = 2 (blue -class (i)). The yellow line indicates the DC conductivity in the absence of heat flows which is also universally bounding σ from below. Although not obvious in the plot, the orange curve does also meet the σ axis at T = 0 with σ(0) ≈ 7.37.
electric conductivity is zero. This is obviously not the case in this system because of the presence of the residual incoherent conductivity discussed above. Here, we operationally define an insulating state as a phase in which dσ/dT > 0. To find the transition points (IMC and MIC), we write (5.35) purely in terms of r 0 , and differentiate the former with respect to the latter. Since the temperature is strictly increasing, we may use the chain rule ∂σ ∂T = ∂σ ∂r 0 to establish a 1-1 relation between each critical point, T k , of σ(T ) and some unique positive real root of ∂σ/∂r 0 , r 0 (T k ) > r e . The three classes are represented in figure 6. As just observed in the context of the DC electric conductivity, these more general gravitational duals open the path for a much richer phenomenology, specially in the thermal transport sector which is strongly modified by the Horndeski terms. In particular, one could think to modify not only the gravity/scalar sectors but also the Maxwell one and to consider more general asymptotics such as Lifshitz geometries, or geometries with hyperscaling violation. A complete analysis of the phenomenology of the dual field theories goes beyond the scope of this work and it will be considered in the near future.

Conclusions
Despite the successes of GR, there are still plenty of reasons to reasonably modify the cornerstone of gravitation, the Einstein-Hilbert action. From the cosmological constant problem and the questionable validity of GR in the very deep infrared (that is, at very large distances), up to quantum gravity corrections and bottom-up holographic applications, looking for robust alternatives is a justified course of action, and in this direction, scalartensor theories are perhaps the most prominent examples that do the job in an absolutely essential manner, i.e., by introducing an additional scalar degree of freedom which may or may not couple non-minimally to gravity. In particular, higher-derivative self-interactions prove to be quite interesting deformations, but they usually involve the risk of ending up with an unstable theory. This endogenous danger can certainly be avoided if the higherderivative contributions somehow cancel out at the level of the field equations, which is exactly the idea behind the theory of Horndeski [29] and its re-emergence in the context of Galileons [68]; by adding appropriate gravitational counterterms, one may suppress the appearance of higher-derivative terms in the equations of motion, keeping them up to second-order in derivatives of the fields.
Still, just putting forth an action principle, although important in its own right, is for various reasons not enough, and one generally wishes to associate a physically meaningful and interesting solution spectrum with it. However, the task of finding exact solutions to a variationally obtained set of field equations is usually a highly non-trivial one, with its complexity depending on the complexity of the model under study and the amount of symmetry being introduced. In this regard, Horndeski theory poses a significant challenge, especially if one is looking for exact black hole solutions. As expected, this road is paved with no-go theorems, dead ends, intractable expressions and naked singularities. Nonetheless, if one manages to overcome these obstacles, a rewarding plethora of highly interesting black hole solutions may be found. Yet, up to today, a large region of the "black" landscape of Horndeski gravity surely remains unexplored, mainly because of the technical difficulties involved, part of which also has to do with the fact that the scalar field is usually taken to depend on the radial coordinate r, a statement compatible with staticity and spherical symmetry.
In the current work, we consider the Einstein-Maxwell theory deformed by the truncated version (2.5) of shift-and reflection-symmetric Horndeski gravity. We show that for backgrounds of the form (1.1), the difficulty of having to integrate the scalar equations of motion is completely lifted, as the latter are identically satisfied. Note, however, that this happens at the cost of reducing the full symmetry of (2.3) down to a residual diagonal subgroup, thus, also partially breaking diffeomorphism invariance in the bulk. For this fairly general class of theories, and starting from a minimal set of requirements found in the beginning of section 3, we provide a systematic way of integrating the metric equations of motion which is schematically given in figure 1. This has the form of a solution-generating algorithm, a step-sequence which we divide it into three separate phases. The algorithm is meant as a means of finding exact solutions respecting these requirements; in particular, this means that the integrals in (3.16), (3.17) and (3.18) must have a closed-form expression, otherwise the whole machinery will not have the desired result. The success of this method relies on an "ad hoc" change of variables, from r toX, with no appearance of the unknown metric functions in the inversion formula (3.3) because of the precise backgrounds (1.1). In terms ofX, the field equations become significantly more tractable, and a straightforward integration can be carried out in a convenient sequence of steps (see the detailed phases analysis in the end of section 3).
Some interesting "universal" results can be also reported at an early stage without specifying the form of G 4 and K. Assuming the solutions describe black hole configurations, we show that one can obtain expressions for the Hawking temperature and Wald entropy density also in terms of the Horndeski functions, and some first restrictions can be placed. Additionally, we argue that if zero-temperature configurations exist, their size will not be affected by G 4 , although the latter contributes to their mass. We further solve the field equations near a hypothetical extremal horizon for an AdS 2 ×R 2 IR fixed point, and display the explicit dependence of the radius of AdS 2 on the Horndeski functions. All extremal bulk solutions we present end up interpolating between this IR geometry and some modeldependent UV counterpart.
Due to the form of the algorithm, we argue that it makes sense to classify solutions based on G 4 . A preliminary exploration of the black hole spectrum is performed by putting our recipe to use, where we immediately (and rapidly) disclose new families of planar black holes in example classes G 4 ∝X and G 4 ∝ (−X) m . In both cases, before fixing the form of k-essence, we do a preliminary asymptotic analysis to filter out choices of K that lead to certain unwanted scenarios. Then, taking into account these restrictions, we choose a K, and derive novel exact black branes with axionic hair. Note that the families we find contain the solutions in [31,36]. Apart from AAdS black branes, we derive a particularly interesting Lifshitz-like family in class G 4 ∝ (−X) m . To the best of our knowledge, this is the first family of exact asymptotically Lifshitz black branes in Horndeski theory with z > 1, where z is a function of m.
Finally, in the context of applied bottom-up holography, we start a preliminary investigation of the transport properties of the dual field theory in the case of an AdS boundary. A detailed derivation of the transport matrix is presented for arbitrary Horndeski functions (also in the presence of a time-dependent heat source), where boundary conditions for the perturbations are discussed in detail. We show, using the example of the DC electric conductivity, that this class of models could give rise to a very rich phenomenology (especially in the strongly modified heat sector) which could be potentially relevant to tackle existing open questions in the condensed matter ballpark.
Needless to say, there are several interesting future directions that can be pursued. For example, an appealing course of action would be to reconstruct the algorithm in higher spacetime dimensions, where a general second-order scalar tensor theory would indeed possess a richer quartic sector, taking into account that, for example, higher terms in the Lanczos-Lovelock series could be also introduced for appropriate d. Additionally, either for scalar-tensor or vector-tensor theories, extensions with higher-order field equations can be also studied; this is the case of, e.g., DHOST theories [69,70], in which, although the theory is higher-order, no Ostrogradsky mode propagates due to the degeneracy of the corresponding kinetic matrix. In the latter class, solutions with a non-constant kinetic term for the scalar field are hard to find [71,72], and in this direction, it would be interesting to see if exporting the setup/method in this work to such cases would bear fruit. If solutions are indeed to be found, then one expects that they will be related to the ones disclosed here by means of a disformal transformation, nevertheless with a completely different causal structure.
Moreover, it would be worth exploring if solutions with a purely time-dependent Galileon background of the form φ ∼ t, reminiscent of the ghost-condensate background [73], could be realized without occurring in evident pathological behaviours; perhaps more interesting in this scenario, would be the holographic interpretation of broken time translations and the possible phenomenological implications. Finally, it would be important to provide a complete characterization of the transport properties of the dual field theories and extract further universal features which are independent of a precise choice of the various Horndeski functions. These questions go beyond the investigation performed in this work and are left to be studied in the future.