Excluding static and spherically symmetric black holes in Einsteinian cubic gravity
with unsuppressed higher-order curvature terms

General Relativity (GR) has been successful in describing the gravitational interaction between submillimeter and solar-system scales [1, 2, 3, 4, 5]. At extremely small distances close to the Planck length, however, it is expected that GR is replaced by a quantum theory of gravity with an ultraviolet completion. The attempt for the construction of a power-counting renormalizable gravitational theory was advocated by Stelle [6] by taking into account quadratic-order curvatures to the Einstein-Hilbert action. In string theory, the low-energy effective action also contains quadratic curvature corrections known as a Gauss-Bonnet (GB) term [7, 8]. Moreover, higher-order curvature terms have played a prominent role in conformal field theory with holography [9, 10, 11, 12].

In gravitational theories where the field equations of motion contain derivatives higher than second order in the metric tensor $g_{\mu\nu}$, the system can be unstable due to the emergence of Ostrogradski instability [13, 14] associated with a Hamiltonian unbounded from below. To avoid such a problem, Lanczos [15] and Lovelock [16] constructed gravitational theories with second-order field equations of motion for general, curved backgrounds. Exploiting polynomial functions of the Riemann curvature tensors in four-dimensional spacetime endowed with four-dimensional diffeomorphism invariance, Lovelock showed that the symmetric and divergence-free tensors $A_{\mu\nu}$, which depend on $g_{\mu\nu}$ and its derivatives up to second order, are expressed by a linear combination of the Einstein tensor $G_{\mu\nu}$ and the metric tensor $g_{\mu\nu}$. In spacetime dimensions higher than four, there is the quadratic-order GB curvature scalar affecting the spacetime dynamics. In four dimensions, the GB term corresponds to an Euler density which does not contribute to the field equations of motion. To extract its effect in four dimensions, the GB term needs to be coupled to other degrees of freedom (DOFs) such as scalar or vector fields [17, 18, 19, 20, 21, 22, 23, 24, 25].

If we consider cubic-order curvature combinations constructed from the Riemann tensors, there is also an additional Euler density in the action corresponding to the surface integral in four dimensions [16]. According to the Lovelock theorem, there are no other nontrivial cubic-order terms keeping the field equations of motion up to second order in general, curved backgrounds. If we consider some specific spacetime, however, it is possible to construct nontrivial cubic-order theory whose graviton spectrum shares a similar property to that in GR. On a maximally symmetric background, which includes the Minkowski and de Sitter spacetimes, there is a unique cubic combination ${\cal P}$ which keeps the structure of linearized perturbation equations of motion in Einstein gravity [26]. In this theory, which is dubbed Einsteinian cubic gravity (ECG), the cubic curvature term can give rise to derivatives higher than second order in general, curved backgrounds, so the spacetime dynamics is generally different from that in GR [27].

On a static and spherically symmetric (SSS) background in a vacuum configuration, ECG admits the existence of a black hole (BH) solution whose temporal and radial metric components (denoted as $f$ and $h$, respectively) coincide with each other [28, 29, 30, 31] (see also Refs. [32, 33, 34, 35]). Since $f$ and $h$ are affected by the cubic curvature term, the background geometry is still different from the Schwarzschild BH solution. One can also construct a cubic gravity theory by imposing the condition $f=h$ (up to a time reparametrization freedom in $f$) on the SSS background [36, 37, 38]. This construction allows the existence of the Lagrangian ${\cal P}$ mentioned above as well as two additional cubic combinations ${\cal C}$ and ${\cal C}^{\prime}$ defined in Ref. [36]. As we will see later in Sec. II, both ${\cal C}$ and ${\cal C}^{\prime}$ vanish for $f=h$ and hence only the Lagrangian ${\cal P}$ contributes to the single differential equation for $f~{}(=h)$.

If we apply ECG to the cosmological dynamics on the Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime, the field equations of motion arising from ${\cal P}$ contain derivatives higher than second order. On the FLRW background, there is a unique combination ${\cal P}-8{\cal C}$ that leads to the second-order Friedmann equation [39, 40]. This theory–dubbed cosmological Einsteinian cubic gravity (CECG)– was applied to the dynamics of inflation and late-time cosmological epochs [39, 40, 41, 42, 43, 44, 45]. On an exact de Sitter background, tensor perturbations with two polarized modes propagate as in GR without pathological behavior. If one considers a spatially homogeneous Bianchi type I manifold close to the isotropic de Sitter spacetime, however, there are three dynamical propagating DOFs associated with linear perturbations in the odd-parity sector. In Ref. [46], it was shown that, in the regime of small anisotropies, such theory possesses at least one ghost mode as well as short-time-scale tachyonic instability. Hence CECG cannot be used to describe a viable geometric inflationary scenario.

Given that there is an instability problem of CECG on the anisotropic cosmological background, we are now interested in the linear stability of SSS vacuum solutions in ECG. For this purpose, we do not deal with ECG as a trivial effective field theory (EFT) where the cubic Lagrangian is always strongly suppressed relative to the Einstein-Hilbert term. We consider odd-parity perturbations according to the Regge-Wheeler formulation [47] without necessarily imposing the condition $f=h$. We show that the Lagrangian ${\cal P}$ in ECG gives rise to three propagating DOFs in the odd-parity sector. We find that there is at least one ghost mode and that the squared propagation speed of one of the dynamical perturbations is negative for large multipoles in the regime where the effective mass of perturbations is below the Planck mass, or the cutoff of the theory $M$. The presence of these ghost and Laplacian instabilities excludes the SSS vacuum solutions in ECG with unsuppressed higher-order curvature terms, including the BH solution with $f=h$.

This paper is organized as follows. In Sec. II, we revisit the SSS BH solution present in ECG. In Sec. III, we derive the second-order action of odd-parity perturbations and show how the ghost and Laplacian instabilities arise for the large frequency and momentum modes. Sec. IV is devoted to conclusions.

General cubic gravity theories consist of a combination of cubic products of the Riemann tensor $R_{\alpha\beta\gamma\delta}$, Ricci tensor $R_{\mu\nu}$, and Ricci scalar $R$. We consider the following eight cubic Lagrangians:

	$\displaystyle\mathcal{L}_{1}={{R_{\alpha}}^{\beta}}{}_{\gamma}{}^{\delta}{{R_{\beta}}^{\mu}}{}_{\delta}{}^{\nu}{{R_{\mu}}^{\alpha}}{}_{\nu}{}^{\gamma}\,,\qquad\mathcal{L}_{2}=R_{\alpha\beta}{}^{\gamma\delta}R_{\gamma\delta}{}^{\mu\nu}R_{\mu\nu}{}^{\alpha\beta}\,,\qquad\mathcal{L}_{3}=R_{\alpha\beta\gamma\delta}R^{\alpha\beta\gamma}{}_{\mu}R^{\delta\mu}\,,\qquad\mathcal{L}_{4}=RR_{\alpha\beta\gamma\delta}R^{\alpha\beta\gamma\delta}\,,$
	$\displaystyle\mathcal{L}_{5}=R_{\alpha\beta\gamma\delta}R^{\alpha\gamma}R^{\beta\delta}\,,\qquad\mathcal{L}_{6}=R_{\alpha}{}^{\beta}R_{\beta}{}^{\gamma}R_{\gamma}{}^{\alpha}\,,\qquad\mathcal{L}_{7}=RR_{\alpha\beta}R^{\alpha\beta}\,,\qquad\mathcal{L}_{8}=R^{3}\,.$		(1)

Taking into account the Einstein-Hilbert term $M_{\rm Pl}^{2}R/2$, where $M_{\rm Pl}$ is the reduced Planck mass, we express the total action as

$${\cal S}=\frac{M_{\rm pl}^{2}}{2}\int{\rm d}^{4}x\sqrt{-g}\left(R+\sum_{i=1}^{8}c_{i}\mathcal{L}_{i}\right)\,,$$

(2)

where $g$ is a determinant of the metric tensor $g_{\mu\nu}$, and $c_{i}$’s are constants.

ECG is constructed to possess a transverse and massless graviton spectrum as in GR on a maximally symmetric background [26]. Requiring also that the relative coefficients of different curvature terms are the same in all dimensions, there is a unique combination ${\cal P}$ which is neither trivial nor topological in four dimensions. This nontrivial cubic interaction corresponds to the choice of coefficients $c_{1}=12$, $c_{2}=1$, $c_{5}=-12$, and $c_{6}=8$, i.e.,

$${\cal P}=12\mathcal{L}_{1}+\mathcal{L}_{2}-12\mathcal{L}_{5}+8\mathcal{L}_{6}\,.$$

(3)

In general curved spacetime including the SSS background as well as the FLRW background, the equations of motion following from the Lagrangian ${\cal P}$ are higher than the second order. In such cases, higher-order derivatives can induce extra DOFs in comparison to those in GR.

In Ref. [36], the authors took a different approach to the construction of cubic gravity theories (dubbed generalized quasi-topological gravity) by demanding that the vacuum SSS solution is fully characterized by a single field equation. In general, the line element on the SSS background is given by

$${\rm d}s^{2}=-f(r){\rm d}t^{2}+h^{-1}(r){\rm d}r^{2}+r^{2}\left({\rm d}\theta^{2}+\sin^{2}\theta\,{\rm d}\varphi^{2}\right)\,,$$

(4)

where $f$ and $h$ are functions of the radial coordinate $r$. If $f(r)$ is proportional to $h(r)$ such that $f(r)=Nh(r)$, where $N$ is a constant, the time and radial components of the field equations coincide with each other. This gives the following constraints among the coefficients in Eq. (2):

$$c_{4}=\frac{3c_{1}-36c_{2}-14c_{3}}{56},\quad c_{5}=-\frac{3c_{1}+48c_{2}+14c_{3}}{7},\quad c_{7}=\frac{6c_{1}+96c_{2}+14c_{3}-21c_{6}}{28},\quad c_{8}=-\frac{3c_{1}+20c_{2}-7c_{6}}{56}.$$

(5)

The six-dimensional Euler density ${\cal X}_{6}$ corresponds to the coefficients $c_{1}=8$, $c_{2}=-4$, $c_{3}=24$, and $c_{6}=-16$, but this is a topological term that does not affect the field equations of motion. One of the remaining three cubic interactions is the Lagrangian ${\cal P}$ given by Eq. (3). There are also two additional terms

	$\displaystyle{\cal C}$	$\displaystyle=$	$\displaystyle\mathcal{L}_{3}-\frac{1}{4}\mathcal{L}_{4}-2\mathcal{L}_{5}+\frac{1}{2}\mathcal{L}_{7}\,,$		(6)
	$\displaystyle{\cal C}^{\prime}$	$\displaystyle=$	$\displaystyle\mathcal{L}_{6}-\frac{3}{4}\mathcal{L}_{7}+\frac{1}{8}\mathcal{L}_{8}\,,$		(7)

which correspond to the coefficients $c_{1}=0$, $c_{2}=0$, $c_{3}=1$, $c_{6}=0$ and $c_{1}=0$, $c_{2}=0$, $c_{3}=0$, $c_{6}=1$, respectively. As we will see below, for $f=h$, the contributions of the Lagrangians ${\cal C}$ and ${\cal C}^{\prime}$ to the field equations of motion vanish. On the FLRW spacetime, the Lagrangian ${\cal P}$ gives rise to the Friedmann equation higher than second-order. The specific combination ${\cal P}-8{\cal C}$ leads to the second-order field equations on the FLRW background [39, 40]. If we apply this latter cubic theory to inflation, the presence of small anisotropies close to the de Sitter background generates the instability of cosmological solutions [46].

In both ECG and CECG, the dynamical DOFs are more than those in GR (two tensor polarizations) around general, curved backgrounds. The property that the propagating DOFs around the maximally symmetric background degenerate to those in GR implies that the disappearing DOFs are, in general, strongly coupled [44, 48]. In such theories, the maximally symmetric background corresponds to a singular surface at which the coefficients of higher-order kinetic terms appearing in the perturbation equations of motion are degenerate. This degeneracy leads to the divergence of couplings of the canonically normalized fields, which is a signal of the strong coupling. On the background different from the maximally symmetric spacetime, some pathological behavior like the instability of perturbations usually arises as it happens for the cosmological background [46, 44]. In this paper, we would like to study whether this is also the case for BHs on the SSS background.

On the SSS background given by the line element (4), lets us consider cubic gravity theories given by the action

$${\cal S}=\frac{M_{\rm pl}^{2}}{2}\int{\rm d}^{4}x\sqrt{-g}\left(R+\alpha{\cal P}+\kappa{\cal C}+\mu{\cal C}^{\prime}\right)\,,$$

(8)

where $\alpha$, $\kappa$, $\mu$ are constants. Then, the quantities (3), (6), and (7) reduce, respectively, to

	$\displaystyle{\cal P}$	$\displaystyle=$	$\displaystyle\frac{1}{f^{4}r^{4}}[3f^{\prime 4}h^{3}r^{2}-3f^{3}h^{\prime}(2f^{\prime\prime}hr+f^{\prime}h^{\prime}r-2f^{\prime}h)(h^{\prime}r-2h+2)-6f^{\prime 2}fh^{2}r\{f^{\prime\prime}hr+f^{\prime}(h^{\prime}r+h-1)\}$		(9)
			$\displaystyle+6rf^{\prime}f^{2}h\{f^{\prime}h^{\prime}(h^{\prime}r+h)+f^{\prime\prime}h(h^{\prime}r+2h-2)\}]\,,$
	$\displaystyle{\cal C}$	$\displaystyle=$	$\displaystyle-\frac{3}{8f^{4}r^{4}}\left(f^{\prime}h-h^{\prime}f\right)^{2}\left[(2f^{\prime\prime}fh-hf^{\prime 2}+f^{\prime}h^{\prime}f)r^{2}-4f^{2}(h-1)\right]\,,$		(10)
	$\displaystyle{\cal C}^{\prime}$	$\displaystyle=$	$\displaystyle\frac{1}{2}{\cal C}\,,$		(11)

where a prime represents the derivative with respect to $r$. If $f$ is equal to $h$, both ${\cal C}$ and ${\cal C}^{\prime}$ vanish. Varying the action (8) with respect to $f$ and $h$, we obtain the third-order differential equations for $h$ and $f$, respectively¹¹1If we relax the assumption of staticity, it can be shown that for this theory the Birkhoff theorem does not hold in general. Therefore, there could be other solutions that in principle have some relevance but will be in general time-dependent. We will not investigate their existence, however, as we will already set strong constraints on the theory just by looking at their static limit.. Setting $f(r)=h(r)$, the two differential equations coincide with each other, giving

$$rf^{\prime}+f-1+\frac{6\alpha}{r^{3}}[rf^{\prime 2}(4f-1)+rf\{r^{2}f^{\prime\prime 2}+4f^{\prime\prime}(f-1)-2rf^{\prime\prime\prime}(f-1)\}+f^{\prime}f(r^{3}f^{\prime\prime\prime}-4r^{2}f^{\prime\prime}-4f+4)]=0\,.$$

(12)

For $\alpha=0$, the solution to Eq. (12) is $f=h=1-r_{h}/r$, where $r_{h}$ is an integration constant corresponding to the horizon radius. For $\alpha\neq 0$, there are corrections to the Schwarzschild metric. At large distances away from the horizon ($r\gg r_{h}$), we derive the solution to Eq. (12) under the following expansion

$$f(r)=h(r)=1-\frac{r_{h}}{r}+\sum_{i=2}c_{i}\left(\frac{r_{h}}{r}\right)^{i}\,,$$

(13)

where $c_{i}$’s are constants. Substituting Eq. (13) into Eq. (12) and computing the coefficients at each order, we obtain $c_{2}=c_{3}=c_{4}=c_{5}=0$, $c_{6}=54\alpha/r_{h}^{4}$, and $c_{7}=-46\alpha/r_{h}^{4}$. Then, the leading-order correction to $f$ arises at sixth order in the expansion (12), so that²²2Here we do not study whether the series converges, as in any case, the solution loses its validity for large enough values of $r$, but we assume that it can be considered as a good approximation in the physical range of $r$ of interest to the real numerical solution having suitable asymptotically flat boundary conditions.

$$f(r)=h(r)=1-\frac{r_{h}}{r}+\frac{54\alpha r_{h}^{2}}{r^{6}}\left[1+{\cal O}\left(\frac{r_{h}}{r}\right)\right]\,.$$

(14)

The solution (14) has been derived by requiring that $f(r)=h(r)$ without explicitly imposing any requirement on the value of $\alpha$. The same solution also follows without assuming the condition $f(r)=h(r)$ in ECG with

$${\cal P}\neq 0\,,\qquad{\cal C}=0\,,\qquad{\cal C}^{\prime}=0\,.$$

(15)

In this case, we write the large-distance solutions as

$$f(r)=1-\frac{r_{h}}{r}+\sum_{i=2}c_{i}\left(\frac{r_{h}}{r}\right)^{i}\,,\qquad h(r)=1-\frac{r_{h}}{r}+\sum_{i=2}d_{i}\left(\frac{r_{h}}{r}\right)^{i}\,.$$

(16)

We substitute Eq. (16) into the third-order differential equations of $f$ and $h$ and obtain the coefficients $c_{i}$ and $d_{i}$. This gives rise to the same coefficients $c_{i}$ derived above, with $c_{i}=d_{i}$, so that the large-distance solution is again given by Eq. (14). Hence the BH solution with $f(r)=h(r)$ generically arises in ECG.

We study the stability of SSS vacuum solutions in ECG given by the action

$${\cal S}=\frac{M_{\rm pl}^{2}}{2}\int{\rm d}^{4}x\sqrt{-g}\left(R+\alpha{\cal P}\right)\,.$$

(17)

In our analysis, we do not restrict ourselves to the EFT regime where $\alpha{\cal P}$ is always suppressed relative to $R$. The modification to the Schwarzschild BH in GR can be significant by allowing for the possibility that $\alpha{\cal P}$ can be as large as $R$. This is analogous to Starobinsky’s model given by the Lagrangian $L=R+\beta R^{2}$ [49], where the cosmic acceleration occurs in the regime $\beta R^{2}\gtrsim R$. If we stick to the EFT regime with $\beta R^{2}\ll R$, one cannot accommodate the physics of inflation driven by the quadratic curvature term.

On the background (4), we decompose the metric tensor into $g_{\mu\nu}=g_{\mu\nu}^{(0)}+h_{\mu\nu}$, where $g_{\mu\nu}^{(0)}$ and $h_{\mu\nu}$ correspond to the background and perturbed parts, respectively. Although we are primarily interested in the stability of the BH solution (14), we do not impose the condition $f=h$ for the background metric from the beginning. Under the rotation in the $(\theta,\varphi)$ plane, the metric perturbations $h_{\mu\nu}$ can be separated into odd- and even-parity modes [47, 50]. Expanding $h_{\mu\nu}$ in terms of the spherical harmonics $Y_{lm}(\theta,\varphi)$, the odd- and even-modes have parities $(-1)^{l+1}$ and $(-1)^{l}$, respectively. In the odd-parity sector, the components of $h_{\mu\nu}$ are given by [51, 52, 53, 54, 55]

	$\displaystyle h_{tt}=h_{tr}=h_{rr}=0\,,\qquad h_{ab}=0\,,$
	$\displaystyle h_{ta}=\sum_{l,m}Q(t,r)E_{ab}\nabla^{b}Y_{lm}(\theta,\varphi)\,,\qquad h_{ra}=\sum_{l,m}W(t,r)E_{ab}\nabla^{b}Y_{lm}(\theta,\varphi)\,,$
	$\displaystyle h_{ab}=\frac{1}{2}\sum_{l,m}U(t,r)\left[E_{a}{}^{c}\nabla_{c}\nabla_{b}Y_{lm}(\theta,\varphi)+E_{b}{}^{c}\nabla_{c}\nabla_{a}Y_{lm}(\theta,\varphi)\right]\,,$		(18)

where $Q$, $W$, and $U$ depend on $t$ and $r$, and the subscripts $a$ and $b$ denote either $\theta$ or $\varphi$. $E_{ab}$ is an antisymmetric tensor with nonvanishing components $E_{\theta\varphi}=-E_{\varphi\theta}=\sin\theta$. In a strict sense, we should write subscripts $l$ and $m$ for the variables $Q$, $W$, and $U$, but we omit them for brevity. Under the gauge transformation $x_{\mu}\to x_{\mu}+\xi_{\mu}$, where $\xi_{t}=0$, $\xi_{r}=0$, and $\xi_{a}=\sum_{l,m}\Lambda(t,r)E_{ab}\nabla^{b}Y_{lm}(\theta,\varphi)$, metric perturbations transform as $Q\to Q+\dot{\Lambda}$, $W\to W+\Lambda^{\prime}-2\Lambda/r$, and $U\to U+2\Lambda$. There are several gauge-invariant combinations like

$$\hat{W}\equiv W-\frac{1}{2}U^{\prime}+\frac{1}{r}U\,,\qquad\hat{Q}\equiv Q-\frac{1}{2}\dot{U}\,.$$

(19)

We fix the residual gauge DOF by choosing

$$U=0\,,$$

(20)

under which $\hat{W}=W$ and $\hat{Q}=Q$.

For the purpose of expanding the action (17) up to second order in odd-parity perturbations, we will focus on the axisymmetric modes ($m=0$) without loss of generality because nonaxisymmetric modes ($m\neq 0$) can be restored under a suitable rotation. For the integrals with respect to $\theta$ and $\varphi$, we exploit the following properties

	$\displaystyle\int_{0}^{2\pi}{\rm d}\varphi\int_{0}^{\pi}{\rm d}\theta\,(Y_{l0,\theta})^{2}\sin\theta=l(l+1)\,,\qquad\int_{0}^{2\pi}{\rm d}\varphi\int_{0}^{\pi}{\rm d}\theta\,\left[\frac{(Y_{l0,\theta})^{2}}{\sin\theta}+(Y_{l0,\theta\theta})^{2}\sin\theta\right]=l^{2}(l+1)^{2}\,,$
	$\displaystyle\int_{0}^{2\pi}{\rm d}\varphi\int_{0}^{\pi}{\rm d}\theta\,\left[\left(\frac{2}{\sin\theta}-\frac{3}{\sin^{3}\theta}\right)(Y_{l0,\theta})^{2}+\frac{3}{\sin\theta}(Y_{l0,\theta\theta})^{2}+\sin\theta\,(Y_{l0,\theta\theta\theta})^{2}\right]=\frac{1}{4}l^{2}(l+1)^{2}(4l^{2}-2l-3)\,.$		(21)

After the integration with respect to $\theta$ and $\varphi$, the second-order action of perturbations, which is expressed in the form ${\cal S}^{(2)}=\int{\rm d}t{\rm d}r\,L$, contains time derivatives such as $\ddot{W}^{2}$, $\dot{Q}^{\prime 2}$, and $\dot{Q}^{2}$. They can be factored out as

$$L_{K}\equiv\frac{3\alpha M_{\rm pl}^{2}\sqrt{h}[2f^{\prime}hr+f(2-2h-rh^{\prime})]l(l+1)}{2f^{5/2}r^{2}}\left(\ddot{W}-{\dot{Q}}^{\prime}+\frac{2\dot{Q}}{r}\right)^{2}\,,$$

(22)

where $L_{K}\in L$. Then, the rest of the Lagrangian $L-L_{K}$ does not contain products like $\ddot{W}\dot{Q}^{\prime}$ and $\ddot{W}\dot{Q}$. We introduce a Lagrange multiplier $\chi$ that helps us to understand the presence and behavior of propagating DOFs. Then, the Lagrangian equivalent to $L$ can be written as $L-L_{K}+\tilde{L}_{K}$ where

$$\tilde{L}_{K}=\frac{3\alpha M_{\rm pl}^{2}\sqrt{h}[2f^{\prime}hr+f(2-2h-rh^{\prime})]l(l+1)}{2f^{5/2}r^{2}}\left[2\chi\left(\ddot{W}-{\dot{Q}}^{\prime}+\frac{2\dot{Q}}{r}\right)-\chi^{2}\right]\,.$$

(23)

Indeed, the variation of (23) with respect to $\chi$ leads to $\chi=\ddot{W}-{\dot{Q}}^{\prime}+2\dot{Q}/r$, so that (23) reduces to (22). We note that the gauge-invariant field $\chi=\ddot{W}-{\dot{Q}}^{\prime}+2\dot{Q}/r$ corresponds to a time derivative of the dynamical perturbation $\tilde{\chi}=\dot{W}-Q^{\prime}+2Q/r$ in GR [47, 53, 55] (up to a free function of $r$). Integrating the Lagrangian $\chi\ddot{W}$ by parts in Eq. (23), we obtain the product $-\dot{W}\dot{\chi}$ between the two first derivatives. As we will see below, there are three propagating DOFs $\vec{{\cal X}}=(W,\chi,Q)$ in this system. In other words, we generally need to give six independent initial conditions to determine the time evolution of $\vec{{\cal X}}$.

After the integration by parts, the total second-order can be expressed in the form ${\cal S}^{(2)}=\int{\rm d}t{\rm d}r\,L$, where

	$\displaystyle L$	$\displaystyle=$	$\displaystyle a_{1}\dot{W}^{2}+a_{2}\dot{Q}^{2}+2a_{3}\dot{W}\dot{\chi}+a_{4}\left(\dot{W}^{\prime}-Q^{\prime\prime}+\frac{2Q^{\prime}}{r}\right)^{2}+a_{5}W^{\prime 2}+a_{6}Q^{\prime 2}+a_{7}W^{2}+a_{8}\chi^{2}+a_{9}Q^{2}$		(24)
			$\displaystyle+a_{10}W^{\prime}\dot{Q}+a_{11}\dot{W}Q^{\prime}+a_{12}\dot{\chi}Q^{\prime}+a_{13}\dot{W}Q+a_{14}\chi\dot{Q}\,,$

with $a_{1},\cdots,a_{14}$ being functions of $r$ alone. Varying the Lagrangian (24) with respect to $W$, $\chi$, and $Q$, we obtain the field equations of motion for these perturbations. They contain the derivatives up to the fourth order, e.g., $\ddot{W}^{\prime\prime}$, $\dot{Q}^{\prime\prime\prime}$.

In the following, we study the propagation of fast oscillating modes by assuming the solutions in the form

$$\vec{\cal X}=\vec{\cal X}_{0}e^{i(\omega t-kr)}\,,\qquad{\rm with}\qquad\vec{\cal X}_{0}=(W_{0},\chi_{0},Q_{0})\,.$$

(25)

Because of staticity, here, the coefficients $\vec{{\cal X}}_{0}$ are supposed to be functions slowly varying in the $r$-direction, so that in the WKB domain they satisfy for instance the condition $|{\vec{{\cal X}}}_{0}^{\prime}|\ll|k\vec{{\cal X}}_{0}|$. For the same reason, we will also suppose that $|\omega^{\prime}|\ll|k\omega|\simeq|\omega^{2}|$. As a consequence, we will also consider the wavenumber $k$ and the multipole $l$ in the ranges $kr_{h}\gg 1$ and $l\gg 1$.

Then, each coefficient in Eq. (24) has the following $l$ dependence:

	$\displaystyle a_{1}=b_{1}l^{4}\,,\quad a_{2}=b_{2}l^{4}\,,\quad a_{3}=b_{3}l^{2}\,,\quad a_{4}=b_{4}l^{2}\,,\quad a_{5}=b_{5}l^{4}\,,\quad a_{6}=b_{6}l^{4}\,,\quad a_{7}=b_{7}l^{6}\,,\quad a_{8}=b_{8}l^{2}\,,$
	$\displaystyle a_{9}=b_{9}l^{6}\,,\quad a_{10}=b_{10}l^{4}\,,\quad a_{11}=b_{11}l^{4}\,,\quad a_{12}=b_{12}l^{2}\,,\quad a_{13}=b_{13}l^{4}\,,\quad a_{14}=b_{14}l^{2}\,,$		(26)

where

	$\displaystyle b_{1}=\frac{3\alpha M_{\rm Pl}^{2}h^{1/2}[fr(2f^{\prime\prime}hr+f^{\prime}h^{\prime}r-6f^{\prime}h)+2f^{2}(h^{\prime}r+2h-2)-f^{\prime 2}hr^{2}]}{4f^{5/2}r^{4}}\,,$
	$\displaystyle b_{2}=\frac{3\alpha M_{\rm Pl}^{2}[(2f^{\prime\prime}f-f^{\prime 2})hr+h^{\prime}f(f^{\prime}r-2f)]}{2f^{7/2}h^{1/2}r^{3}}\,,$
	$\displaystyle b_{3}=\frac{r^{2}}{f}b_{1}-\frac{hr^{2}}{2}b_{2}+\frac{b_{0}}{2fh}\,,\qquad b_{4}=fhb_{3}+\frac{3}{2}b_{0}\,,\qquad b_{5}=f^{2}h^{2}b_{2}-\frac{f}{r^{2}}b_{0}\,,\qquad b_{6}=b_{1}+\frac{2}{hr^{2}}b_{0}\,,$
	$\displaystyle b_{7}=-\frac{b_{5}}{2hr^{2}}\,,\qquad b_{8}=b_{3}\,,\qquad b_{9}=\frac{f}{2r^{2}}b_{2}\,,\qquad b_{10}=-\frac{3\alpha M_{\rm Pl}^{2}\sqrt{h}(3f^{\prime}h+h^{\prime}f)}{f^{3/2}r^{3}}\,,$
	$\displaystyle b_{11}=-2b_{1}-2fhb_{2}-b_{10}-\frac{b_{0}}{hr^{2}}\,,\qquad b_{12}=-2b_{3}\,,$		(27)

with

$$b_{0}=\frac{3\alpha M_{\rm Pl}^{2}h^{3/2}(f^{\prime}h-h^{\prime}f)}{f^{3/2}r}\,.$$

(28)

We do not show the explicit forms of $b_{13}$ and $b_{14}$, as they are not needed in the following discussion.

Exploiting the approximation of large values of $\omega$, $k$, and $l$ in the field equations of motion, we ignore terms proportional to $i\omega$, $ik$ and $l$ relative to those proportional to $\omega^{2}$, $k^{2}$ and $l^{2}$. Then, the perturbation equations of motion following from the Lagrangian (24) can be expressed in the form

$${\bm{A}}\vec{{\cal X}}_{0}^{\rm T}=0\,,$$

(29)

where ${\bm{A}}$ is the $3\times 3$ matrix whose components are given by

	$\displaystyle A_{11}=2l^{2}\left[(b_{4}k^{2}+b_{1}l^{2})\omega^{2}+b_{5}k^{2}l^{2}+b_{7}l^{4}\right]\,,\qquad A_{22}=2l^{2}b_{8}\,,\qquad A_{33}=2l^{2}\left(b_{2}l^{2}\omega^{2}+b_{4}k^{4}+b_{6}k^{2}l^{2}+b_{9}l^{4}\right)\,,$
	$\displaystyle A_{12}=A_{21}=2l^{2}b_{3}\omega^{2}\,,\qquad A_{13}=A_{31}=l^{2}(2b_{4}k^{3}-b_{10}kl^{2}-b_{11}kl^{2})\omega\,,\qquad A_{23}=A_{32}=-l^{2}b_{12}k\,\omega\,.$		(30)