Gravitational Decoupling in Higher Order Theories

: Gravitational decoupling via the Minimal Geometric Deformation (MGD) approach has been used extensively in General Relativity (GR), mainly as a simple method for generating exact anisotropic solutions from perfect ﬂuid seed solutions. Recently this method has also been used to generate exact spherically symmetric solutions of the Einstein-scalar system from the Schwarzschild vacuum metric. This was then used to investigate the effect of scalar ﬁelds on the Schwarzschild black hole solution. We show that this method can be extended to higher order theories. In particular, we consider fourth order Einstein–Weyl gravity, and in this case by using the Schwarzschild metric as a seed solution to the associated vacuum ﬁeld equations, we apply the MGD method to generate a solution to the Einstein–Weyl scalar theory representing a hairy black hole solution. This solution is expressed in terms of a series using the Homotopy Analysis Method (HAM).


Introduction
Obtaining exact physically relevant solutions of Einstein's field equations is by no means a straightforward task. The spherically symmetric perfect fluid analytic solutions of Einstein's equations have been derived and analyzed in detail [1,2]. However, as soon as one relaxes the assumption of a perfect fluid source and considers more realistic sources, then the problem of finding analytic solutions becomes intractable. Fortunately there are a number of solution generating techniques in General Relativity (GR) [3][4][5], where one starts with a seed metric, which is often a vacuum or a perfect fluid solution, and use these to generate other solutions with more complex forms of the energy momentum tensor that are more applicable to realistic scenarios. Recently, the Minimal Geometric Deformation (MGD) technique has been proposed as a novel and simple approach to decoupling gravitational sources in General Relativity, which could lead to new static and spherically symmetric solutions having sources that are more realistic than the ideal perfect fluid. This method was originally introduced in the context of the Randall-Sundrum braneworld [6,7] and has been used to generate brane-world configurations [8,9] from general relativistic perfect fluid seed solutions. The method was later extended to investigate new black hole solutions [10,11] and to generate physically meaningful interior stellar solutions [12] (for some other more recent applications see Refs. [13][14][15][16][17][18][19][20][21]). Gravitational decoupling has also been studied in modified theories of gravity such as f (R, T) [22,23], Gauss-Bonnet gravity [24], f (G) gravity [25] and Rastall gravity [26].
The main feature of this method is that one can start with a known seed solution of Einstein's field equations having a relatively simple sourceT µν , such as, for example, a perfect fluid solution, and add to this a second source θ µν such that the total energy momentum tensor is T µν =T µν + θ µν , where θ µν may be a generic anisotropic source. Through the MGD the new Einstein field equations for T µν can then be separated into two sets of equations, that is, the original system of Einstein's field equations forT µν leading to the seed solution, and another system of equations for the additional source θ µν . The solution of the second system of equations leads to the complete solution corresponding to the total source T µν . This process can then be repeated by adding additional sources to T µν to finally end up with an intricate source that describes a realistic scenario. The converse can also be applied. In this case, one starts with a complicated energy momentum tensor which can be split into simpler components. The corresponding Einstein field equations for the whole system can be reduced to separate systems of equations that can be solved separately. In both cases, the method works provided that there is no interaction between the individual sources, that is, the energy momentum tensors of all the sources are conserved separately, such that ∇ νT µν = ∇ ν θ µν = 0. Moreover, the gravitational decoupling is attained due to the time independence and spherical symmetry of the system. The MGD method has also been used to investigate solutions of Einstein-scalar gravity [27]. In this case, the additional source θ µν corresponds to the energy momentum tensor of a minimally coupled scalar field. It has been shown [28] that, for some perfect fluid configurations, adding a scalar field will make the complete system unstable.
The main objective of this study is to extend the MGD technique to higher derivative theories of gravity and in particular to pure Einstein-Weyl gravity. It is a known fact that Einstein's general relativity is non-renormalisable. In higher derivative theories of gravity, the Einstein-Hilbert action is changed by including quadratic curvature invariants [29] in the action so as to get a renormalisable theory. The most general action in this case takes the form: where g is the determinant of the metric tensor g µν , β, δ and γ are constants, R is the Ricci scalar and C µνρσ is the Weyl tensor. In general when including higher order terms in the action, one has to consider the possibility of ghost-like modes in the theory. However these do not always render the theory irrelevant, as shown in Ref. [30]. Therefore these higher derivative extensions of general relativity can be considered as alternative and effective theories of gravity. Besides the fact that these arise naturally in the string theory approach to gravity, they are mainly used in cosmology to generate geometric dark energy models without the need to use a cosmological constant, unlike the Λ Cold Dark Matter (ΛCDM) model of General Relativity. Any vacuum solution of Einstein's general relativity is also a solution to the higher derivative theory in (1), but the converse is not necessarily true. In fact, the higher derivative theory admits a richer set of vacuum solutions than General Relativity. For example, by taking δ = 0 in the action (1), which is referred to as pure Einstein-Weyl gravity, the authors in Refs. [31,32] obtained a static and spherically symmetric asymptotically flat vacuum black hole solution that is different than the Schwarzschild solution. This solution was obtained numerically. The authors also looked at its thermodynamical properties and showed that it satisfies the first law, and for a given mass its entropy is always less than the entropy of the Schwarzschild black hole with the same mass. It has been shown (see Ref. [33]) that any static black hole solution of (1) must have vanishing Ricci scalar, and so in this case the term quadratic in R in the action (1) makes no contribution to the field equations.In this paper we consider the case of pure Einstein-Weyl gravity and include a source θ µν corresponding to a minimally coupled scalar field. Due to the lack of non-trivial analytical solutions in Einstein-Weyl gravity we choose the Schwarzschild metric as our seed solution, since as mentioned earlier, all vacuum solutions of GR are also solutions of Einstein-Weyl gravity. Using the MGD method we obtain the field equations for θ µν , and these are solved using the Homotopy Analysis Method (HAM), which results in an approximate analytic solution representing a hairy black hole. The Homotopy Analysis Method (HAM) [34] is a very useful method for obtaining analytical approximate solutions to various nonlinear differential equations, including also systems of nonlinear equations, arising in many different areas of science and engineering. Finally, the effect of the scalar field on the local black hole geometry is analyzed. The paper is organized as follows. In Section 2 we review the gravitational decoupling by the MGD method in General Relativity. Then in Section 3 we extend this method to pure Einstein-Weyl gravity and obtain the field equations for the Einstein-Weylscalar system. The Homotopy Analysis Method (HAM) is briefly introduced in Section 4, and then in Section 5 this is used to obtain an analytic approximate solution to the system of equations in Einstein-Weyl-scalar gravity. The results are summarized and discussed in the Conclusion. In this paper we will use units such that 8πG = c = 1.

Minimal Geometric Deformation in General Relativity
Consider Einstein's field equations for an interior stellar distribution with a source T µν , which consists of a perfect fluid and an additional source θ µν , such that where ρ is the fluid density, p is the isotropic fluid pressure, u µ is the four-velocity of the fluid and α is a coupling constant. The additional source θ µν may contain new scalar, vector or tensor fields, and will in general produce the anisotropies in the self-gravitating system. The total energy momentum tensor T µν is conserved and, since we assume that there is no exchange of energy between the perfect fluid and the additional source represented by θ µν , we can write Then, for the static spherically symmetric metric in Schwarzschild coordinates, the Einstein's field equations are given by where indicates a derivative with respect to the areal coordinate radius r. The conservation equation ∇ ν T µν = 0 takes the form: h ρr + h rp + 2hrp − αh rθ t t + αh rθ r r + 2αhr(θ r r ) + 4αhθ r r − 4αhθ θ θ = 0.
As can be seen from the field equations in (7), the additional source θ µν introduces an anisotropy inside the stellar distribution such that the effective energy density is ρ = ρ − αθ t t , while the effective radial and tangential pressures are given byp r = p + αθ r r andp t = p + αθ θ θ , respectively. So the field equations in (7) contain five unknown parameters, namelyρ,p r ,p t and the metric functions h(r), f (r), and therefore the system is under determined.
The geometric deformation of the radial component of the metric f (r), given by [28,35] where µ(r) is the deformation function, leads to a gravitational decoupling of the Einstein's field Equation (7) into two separate systems of equations, namely: and Similarly, the conservation equation ∇ µ T µν = 0 leads to One can realize that the systems in (10) are the field equations for a perfect fluid solution (ρ, p, h(r), f (r)) (referred to as the seed solution), and (11) is the corresponding conservation equation. So given a seed solution, one can generate an anisotropic solution by solving the system of Equation (11) together with the conservation Equation (12). Since (12) follows from (11), the system of Equation (11) is still under determined, and therefore to obtain a solution one would need to impose a constraint, such as θ r r (r) = p(r). This would lead to the anisotropic solution (ρ(r),p r (r),p t (r)) with metric functions (h(r), f (r) + αµ(r)).
If the additional source is chosen to be a minimally coupled scalar field with scalar potential V(φ), whose action is given by then the corresponding energy momentum tensor is given by The conservation Equation (12) becomes the Klein-Gordon equation where ≡ ∇ µ ∇ µ . In this case the system of Equation (11) are given by which can be solved directly for the unknown functions (µ, φ, V). In other words, when using the MGD method to generate Einstein-scalar solutions, one does not need to impose an additional constraint such as θ r r (r) = p(r), which is required for anisotropic fluid solutions. This method has been used to investigate the effect of a scalar field on the Schwarzschild vacuum solution [27].

Einstein-Weyl Gravity
As remarked earlier, the action of pure Einstein-Weyl gravity in a vacuum is given by (1) with δ = 0. Using units in which γ = 1, the equations of motion obtained by varying the action with respect to the metric g µν , are given by where the Bach tensor B µν is expressed in terms of the Weyl tensor by ;ν = 0), and so taking the trace of (17) gives R = 0. The Schwarzschild metric is also a solution of Einstein-Weyl gravity. However, as shown in Ref. [31], this is not the only static and spherically symmetric vacuum solution of Einstein's Weyl gravity.
For the static and spherically symmetric metric (5) the vacuum field equations are given by: and Due to the inclusion of the Weyl squared term in the action, one would expect that the corresponding field Equation (17) for the spherically symmetric metric (5) contains third and fourth order derivatives of the metric functions f (r) and h(r). However, as was done in this case, such higher order derivatives can be eliminated by taking a combination of the field equations such that the resulting Equations (18) and (19) now contain only second order derivatives. Yet they are still too complicated to be solved analytically unlike the corresponding field equations in general relativity. Therefore, in such cases, one normally reverts to numerical methods. This was was done in Ref. [31] (see also [36]), where the authors obtained a Schwarzschild-like asymptotically flat numerical black hole solution. Fixing the event horizon at r = r 0 , the metric functions h(r) and f (r) can be represented by the following Taylor expansions close to the horizon: where c is an arbitrary scaling factor for the time coordinate, chosen such that it matches with the proper time for a remote observer, that is, lim r→∞ h(r) = 1. Then, substituting (20) in the field Equations (18) and (19), one can express the coefficients h i , f i i ≥ 2 in terms of the free parameters r 0 and f 1 , so that and and similarly for the remaining coefficients. Fixing the horizon radius r 0 and using the above expansions (20) to obtain the initial data at r i (where r i is taken to be very close to the horizon radius r 0 ), the field Equations (18) and (19) can be numerically integrated by starting the integration at the distance r i . The shooting method is used in order to find the appropriate values of the parameters c and f 1 in (20) such that the solution is asymptotically flat up to a large distance of r = r f ∼ 40r i . The authors of Ref. [31] showed that there is a certain value of the horizon radius r min 0 which depends on the parameter β, such that for r 0 > r min 0 , the theory admits precisely one static black hole solution besides the Schwarzschild solution. So, for example, letting β = 1/2 they found that r min 0 ≈ 0.876. Figure 1 shows such a non-Schwarzschild black hole for β = 1/2 and r 0 = 1.

Gravitational Decoupling in Einstein-Weyl Gravity
To obtain gravitational decoupling of the field equations in Einstein-Weyl gravity we use the same approach as in GR. We take a minimally coupled scalar field φ as the source, such that the action for Einstein-Weyl-scalar gravity is given by The corresponding field equations are and the conservation equation for the field φ is identical to that used in GR, namely where φ = g µν φ ;µν . Starting with the general static and spherically symmetric metric (5) we apply the minimal geometric deformation and substitute the deformed metric in the above field Equation (24). As in the GR case, we find that this would lead to a decoupling of the field equations into two systems. The first set contains the vacuum field equations of Einstein-Weyl gravity for the metric functions f (r) and h(r) given above by (18) and (19), which were solved numerically in Ref. [31] to obtain the non-Schwarzschild black hole solution shown in Figure 1. The other set is a system of differential equations for the deformation function µ(r), the scalar field φ(r) and its potential V(φ). These equations are also second order and are given by and The above system of equations can be solved numerically for φ, µ, V using the non-Schwarzschild vacuum solution ( f (r), h(r)) obtained in Ref. [31] and displayed in Figure 1. This will generate a numerical hairy black hole solution to the Einstein-Weyl-scalar theory. In our case we decide to use the simple Schwarzschild metric as our seed solution instead. As explained earlier, this metric is also a solution of the vacuum field Equations (18) and (19) in Einstein-Weyl gravity. Then, by solving the above system of equations we can obtain a hairy black hole solution, which can be represented analytically using the Homotopy Analysis Method. This method is explained in detail in the next section.

Homotopy Analysis Method
In this section, we give a brief overview of the basics of the Homotopy Analysis Method (HAM). Finding exact solutions of non-linear differential equations is not always possible. There are various techniques [37] that have been developed to obtain analytical approximate solutions to non-linear differential equations or systems of non-linear differential equations. Among these, one can mention the power series method, perturbation techniques, variational iteration method and sinc-collocation method, just to name a few. Another method that has been utilized extensively during the last two decades for various mathematical and physical problems is the Homotopy Analysis Method (HAM), which was first introduced by Liao in 1992 [34,[38][39][40][41][42] (for various applications of the HAM see [43][44][45][46][47][48][49][50][51][52][53][54][55]). Most of these approximation techniques will generate a series solution to a non-linear differential equation or a system of non-linear differential equations. In this regard, HAM has the unique property that the region and the rate of convergence of the series solution is not just dependent on the choice of the initial approximation, auxiliary linear operator and an auxiliary function, but it can be effectively controlled using a convergence control parameter. This added freedom comes with its own problems, due to the fact that currently there are no mathematical theorems which give the values of these parameters for a given problem. Despite this, there are a number of guidelines (see for example [38,56,57]) that are useful when choosing the right parameters in fairly general situations.
Despite its popularity in many areas of science and engineering, the application of the HAM in General Relativity and gravitation in general has been very limited (see [58] for a recent application). In General Relativity and particularly in higher order theories of gravity such as Einstein-Weyl gravity, the field equations are highly nonlinear and in most cases do not lead to an exact solution for the metric tensor g µν , which means that, in such cases, the only option is to revert to numerical methods. However, a numerical solution can be achieved given a set of initial conditions (or boundary conditions in the case of a boundary value problem) and for fixed values of the parameters that exist in the field equations. So although this can be used for further numerical analysis of the given problem, one would not be able to obtain a clear picture of the dependence of the metric on the parameters used in the system. Therefore, in the absence of an exact solution, it would be desirable to get at least an analytic approximate solution.
To explain the basics of the HAM consider a system of n nonlinear differential equations which we will write in the compact form, where N i are nonlinear operators, y i (t) are the unknown functions and t is the independent variable. For this system we consider the homotopy [38], where L is an auxiliary linear operator, q ∈ [0, 1] is an embedding parameter, H i are auxiliary functions, h i = 0 are the convergence control parameters, y i0 (t) are the initial guesses of y i (t) and φ i (t; q) are unknown functions. Without loss of generality (see for example [56]) one may set the auxiliary functions H i (t) = 1. As mentioned earlier, one has the freedom to assign the auxiliary linear operator, initial approximations and the convergence control parameters. These choices affect the convergence of the solutions and therefore the computational efficiency of the method. As mentioned already, the ways in which these can be chosen are discussed in Refs. [38,56,57]. It is evident from (31) that as q increases from 0 to 1, the solution φ i (t; q) varies from the initial guess y i0 (t) to the actual solution y i (t). For this reason Equation (31) serves as the zeroth order deformation equation. Expanding φ i (t; q) as a Taylor series with respect to q, gives where The choice of the initial approximation, linear operator, and convergence control parameter should lead to convergence of the above Taylor expansion at q = 1 such that One way to choose the parameters h i is to plot the so-called h-curves which represent the graphs of the dependent variables y i (t) or their derivatives calculated at specific values of the independent variable t (normally taken to be t = 0), against h i . If (34) is convergent, a horizontal line segment will be obtained in the corresponding h-curve. Moreover, if a value of h i is chosen from this region, it can be shown that (34) converges (see Ref. [38] for details). The initial approximations y i0 (t) are usually taken to be the simplest functions satisfying the boundary conditions of the system of differential equations [57] and the linear operators are normally taken to be the highest order derivatives of the dependent variables in the differential equations [56].
The functions y im (t) in (34) are obtained by solving the so called mth-order deformation equation. This is obtained by differentiating (31) m times with respect to the parameter q and then letting q = 0 and finally dividing by m!. This is given by where and The boundary conditions for the above differential equation are obtained from the original problem and since it is linear it can easily be solved for y im (t). If we define the partial sum y M i (t) by then y M i (t) will serve as the M-th order approximation to the solution (34).

Hairy Black Hole in Einstein-Weyl Gravity
Substituting the expression for the potential in (29) into (27) and (28) and changing to the variable x = 1 − r 0 /r gives the following second order nonlinear differential equations, and where now indicates a derivative with respect to x. The initial approximations are taken as where Q is a constant representing the scalar charge on the event horizon x = 0 (r = r 0 ). The chosen approximations satisfy the boundary conditions for the above system of differential equations, in the sense that µ 0 vanishes on the horizon and at spatial infinity where x = 1, and φ 0 also vanishes asymptotically. Since both differential equations are second order, the auxiliary linear operators (with the variable t now replaced by x) are given by L(φ i (x; q)) = ∂ 2 φ i (x; q)/∂x 2 , i = 1, 2. The non-linear operators N i (φ i (x; q)) are obtained from the nonlinear operators used in (39) and (40) after expressing these in a form such that each differential equation contains second derivatives of either µ or φ, but not both. As already mentioned in the previous section, without loss of generality we may set the auxiliary functions H i (x) = 1. The convergence control parameters h i are obtained from the h-curves as explained in the previous section. Substituting the initial approximations (41) in the deformation Equation (35) and solving recursively using the boundary conditions y im (0) = y im (1) = 0, we get the series solution (38) for the functions µ(x) and φ(x). In our case we considered third order approximations (M = 3 in (38)) for both functions. So, letting α = 1, β = 1/2 and r 0 = 1, Q = 1 the h-curves µ (0) vs h and φ (0) vs h are shown in Figure 2a,b. From the two h-curves, it is clear that both flat sections are centred at around h = 5/8. Hence, using this value for the convergence control parameter for both differential equations, we can write down the series approximations for the scalar field and deformation function. Due to the complexity of these expressions, these are included in Appendix A, where both functions are expressed again in terms of the radial coordinate r. Letting r 0 = 1 for simplicity, these can be plotted as shown in Figure 3a  Since we take Q = 1, we see that φ(r 0 ) = 1. Both functions vanish for large r so that the generated solution remains asymptotically flat like the seed Schwarzschild metric. The metric function g rr = f (r) + αµ(r) for the generated metric is shown in Figure 4, along with the corresponding Schwarzschild metric function f (r) = 1 − r 0 /r, represented by the dotted curve. Note that the metric function h(r) = 1 − r 0 /r will be unchanged by the MGD method, and so this is common for the seed Schwarzschild metric and the generated solution. The potential function V(r) can be obtained by substituting the expressions for µ(r) and φ(r) in (29) along with the metric functions f (r) = h(r) = 1 − r 0 /r for the Schwarzschild seed metric. This is shown in Figure 5. As expected, the potential vanishes at infinity. Finally, in Figure 6 we show the Kretschmann scalar κ = R αβµγ αβµγ , which quantifies the spacetime curvature outside the black hole, for the generated solution compared with the Schwarzschild case shown by the dotted curve.

Discussion and Conclusions
Gravitational decoupling by the minimal geometric deformation method is a simple yet powerful technique for generating anisotropic solutions from simple perfect fluid seed metrics. This has been applied extensively in GR to obtain physically acceptable solutions describing the interior geometry of stellar distributions. It has some limitations though. It can only be applied to spherically symmetric systems, and if there are multiple sources, then these have to be decoupled such that no energy exchange takes place between them. Moreover to be able to solve the field equations to obtain the deformation function µ(r), one needs to adopt an additional constraint such as θ r r (r) = p(r). This constraint is not required when the additional source θ µν represents a scalar field, and therefore the method is particularly suited for obtaining scalar field extensions of known seed metrics. This would be particularly interesting for studying hairy black holes, that is, black holes surrounded by scalar hair, as has been done recently in Ref. [59] in GR.
The MGD method is not restricted to General Relativity and in fact it has already been extended to a few modified theories of gravity. In this paper we extend the method to Einstein-Weyl gravity and we show that, although here the field equations are much more complex and in general contain fourth order derivatives of the metric functions, one can still achieve gravitational decoupling by following the same steps as in GR. Obviously the same restrictions apply, that is, the spacetime should be static and spherically symmetric and no energy exchange between multiple sources is allowed. Due to the complexity of the field equations, it is very difficult to find non-trivial exact perfect fluid or even vacuum solutions of Einstein-Weyl gravity [60]. So to demonstrate the application of this method we use the Schwarzschild metric as our seed solution of Einstein-Weyl gravity and add a minimally coupled scalar field source. The HAM is then applied to the resulting field equations to obtain analytic approximations for the deformation function µ(r) and scalar field φ(r). The boundary conditions are chosen such that the generated scalar black hole solution is asymptotically flat as shown in Figures 3 and 4. The third order approximations for µ(r) and φ(r) are listed in Appendix A.
There are other modified theories of gravity where the MGD method may be applied to generate scalar field solutions as described in this paper. Another case that has not yet been considered in the literature is massive gravity, which attributes a mass to the graviton field. An example of such a theory is that developed by de Rham, Gabadadze, and Tolley [61,62], commonly referred as dRGT massive gravity. This is particularly interesting because the field equations in dRGT gravity do not include higher derivative terms and so the theory is ghost-free. Its action is given by where κ 2 = 8πG, R is the Ricci scalar, U is the self-interacting potential for gravitons having mass m g , φ α are the Stückelberg scalar fields and as usual S m is the part of the action corresponding to the matter content. It has been shown [63] that this theory has an exact static and spherically solution of the form in (5) with where the constants Λ, γ and η are expressed in terms of the graviton mass m g , such that when m g = 0 the theory reduces to GR and the above solution becomes the Schwarzschild black hole metric. When γ = η = 0, the metric becomes the Schwarzschild-de Sitter black hole solution, which has been shown to be stable [64] in this theory. The cosmological-like term is useful for explaining the current accelerated expansion of the universe [65]. Unlike other theories, like GR for example, this term arises from the graviton mass itself and is not included by hand in the action of the theory. Moreover in this case, apart from the fact that the massive graviton behaves like a kind of dark energy on the large scale, it can also mimic the dark matter halo on the galactic scale through the linear term γr in the above metric, which can be used to fit the rotational curves of most galaxies [66]. One should mention [67] that this metric is also a solution in other theories, such as conformal Weyl gravity [68], f (R) gravity [69] and general relativity [70], and although in each case the metric represents the same manifold at a geometric level, the physical meaning of the individual parameters m, Λ, γ and η is different. One can try to apply the MGD method to massive gravity in the same way that it has been applied before to GR and to Einstein-Weyl gravity as shown in this paper. So, starting with the general static and spherically symmetric metric in (5), one can apply the minimal geometric deformation f (r) → f (r) + αµ(r) and check whether the field equations corresponding to (42) with S m = 0 decouple. If gravitational decoupling of the field equations can be achieved in dRGT theory, one can use the MGD method to generate new scalar solutions in this theory by using (43) as the seed solution and letting S m = S φ in the action (42), where S φ is the action for a minimally coupled scalar field given by (13), and φ is assumed to be conserved independently, that is, it satisfies Equation (25). This will be explored in a future publication.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: MGD Minimal Geometric Deformation HAM Homotopy Analysis Method GR General Relativity