Instability of scalarized compact objects in Einstein-scalar-Gauss-Bonnet theories

We consider the Horndeski theories Horndeski (1974); Deffayet et al. (2009b); Kobayashi et al. (2011) whose action is composed of the four-independent parts

$\displaystyle S$

$\displaystyle=$

$\displaystyle\int d^{4}x\sqrt{-g}{\cal L}=\int d^{4}x\sqrt{-g}\sum_{i=2}^{5}{\cal L}_{i},$

(1)

with the Lagrangian densities given by

	$\displaystyle{\cal L}_{2}$	$\displaystyle:=$	$\displaystyle G_{2}(\phi,X),$
	$\displaystyle{\cal L}_{3}$	$\displaystyle:=$	$\displaystyle-G_{3}(\phi,X)\Box\phi,$
	$\displaystyle{\cal L}_{4}$	$\displaystyle:=$	$\displaystyle G_{4}(\phi,X)R+G_{4X}(\phi,X)\left[\big{(}\Box\phi\big{)}^{2}-\left(\phi^{\alpha\beta}\phi_{\alpha\beta}\right)\right],$
	$\displaystyle{\cal L}_{5}$	$\displaystyle:=$	$\displaystyle G_{5}(\phi,X)G_{\mu\nu}\phi^{\mu\nu}-\frac{1}{6}G_{5X}(\phi,X)\left[(\Box\phi)^{3}-3\Box\phi\left(\phi^{\alpha\beta}\phi_{\alpha\beta}\right)+2\phi_{\alpha}{}^{\beta}\phi_{\beta}{}^{\rho}\phi_{\rho}{}^{\alpha}\right],$		(2)

where $g_{\mu\nu}$ is the spacetime metric, $R$ and $G_{\mu\nu}$ are the Ricci scalar and Einstein tensor associated with $g_{\mu\nu}$, respectively, $\phi$ is the scalar field, $\phi_{\mu}=\nabla_{\mu}\phi$, $\phi_{\mu\nu}=\nabla_{\mu}\nabla_{\nu}\phi$, and so on, with $\nabla_{\mu}$ being the covariant derivative associated with the metric $g_{\mu\nu}$, $X$ represents the canonical kinetic term $X:=-(1/2)g^{\mu\nu}\phi_{\mu}\phi_{\nu}$, and $G_{i=2,3,4,5}(\phi,X)$ are free functions of $\phi$ and $X$.

The Einstein-scalar-GB theory

$\displaystyle{\cal L}_{\rm sGB}=\frac{1}{2\kappa^{2}}\left[R+\zeta X+\xi(\phi)R_{\rm GB}^{2}\right],$

(3)

where $\kappa^{2}(=8\pi G)$ denotes the gravitational constant, $\zeta$ is a constant, and $\xi(\phi)$ is the coupling function of the scalar field $\phi$ to the GB invariant

$\displaystyle R_{\rm GB}^{2}:=R^{2}-4R^{\alpha\beta}R_{\alpha\beta}+R^{\alpha\beta\mu\nu}R_{\alpha\beta\mu\nu},$

(4)

is known as subclass of the Horndeski theories (1) with the following choice of the coupling functions Kobayashi et al. (2011); Langlois et al. (2022)

	$\displaystyle G_{2}(\phi,X)$	$\displaystyle=$	$\displaystyle\frac{1}{2\kappa^{2}}\zeta X+\frac{4}{\kappa^{2}}\xi^{(4)}(\phi)X^{2}\left(3-\ln|X|\right),$		(5)
	$\displaystyle G_{3}(\phi,X)$	$\displaystyle=$	$\displaystyle\frac{2}{\kappa^{2}}\xi^{(3)}(\phi)X\left(7-3\ln|X|\right),$		(6)
	$\displaystyle G_{4}(\phi,X)$	$\displaystyle=$	$\displaystyle\frac{1}{2\kappa^{2}}+\frac{2\xi^{(2)}(\phi)}{\kappa^{2}}X(2-\ln|X|),$		(7)
	$\displaystyle G_{5}(\phi,X)$	$\displaystyle=$	$\displaystyle-\frac{2}{\kappa^{2}}\xi^{(1)}(\phi)\ln|X|,$		(8)

where $\xi^{(n)}:=\partial^{n}\xi(\phi)/\partial\phi^{n}$. We choose $\zeta>0$ so that the scalar field has the correct sign of the kinetic term if the GB coupling vanishes.

We assume the static and spherically symmetric spacetime

$\displaystyle ds^{2}$

$\displaystyle=$

$\displaystyle-f(r)dt^{2}+\frac{dr^{2}}{h(r)}+r^{2}\gamma_{ab}d\theta^{a}d\theta^{b}.$

(9)

and the static scalar field

$\displaystyle\phi=\phi(r),$

(10)

where $f(r)$, $h(r)$, and $\phi(r)$ are the pure functions of the radial coordinate $r$, $\gamma_{ab}$ represents the metric of the unit two-sphere, and the coordinates $\theta^{a}$ run over the angular directions. Substituting the ansatz (9) and (10) into the Lagrangian (3) and varying it with respect to $f$ and $h$, we obtain the equations of motion for $f$ and $h$ which are given by

		$\displaystyle 2\left(r+2(1-3h)\xi^{(1)}\phi^{\prime}\right)h^{\prime}+\frac{1}{2}\left(-4+h(4+\zeta r^{2}\phi^{\prime 2})\right)-8h(h-1)\left(\phi^{\prime 2}\xi^{(2)}+\phi^{\prime\prime}\xi^{(1)}\right)=0,$		(11)
		$\displaystyle 2\left(r+2(1-3h)\xi^{(1)}\phi^{\prime}\right)f^{\prime}-\frac{f}{2h}\left(4+h(-4+\zeta r^{2}\phi^{\prime 2})\right)=0.$		(12)

Similarly, varying the Lagrangian (3) with respect to the scalar field $\phi$, we obtain the scalar field equation of motion as

$\displaystyle\phi^{\prime\prime}+\frac{1}{2}\left(\frac{4}{r}+\frac{f^{\prime}}{f}+\frac{h^{\prime}}{h}\right)\phi^{\prime}+\frac{4}{\zeta r^{2}f}(-1+h)\xi^{(1)}f^{\prime\prime}(r)-\frac{2f^{\prime}\xi^{(1)}}{\zeta r^{2}f^{2}h}\left(h^{2}f^{\prime}+fh^{\prime}-f^{\prime}h-3fhh^{\prime}\right)=0.$

(13)

Solving Eq. (12) with respect to $h^{-1}$, we obtain

$\displaystyle\frac{1}{h}=\frac{1}{8f}\left(4f^{\prime}(r+2\xi^{(1)}\phi^{\prime})+f(4-\zeta r^{2}\phi^{\prime 2})\pm\sqrt{-384ff^{\prime}\xi^{(1)}\phi^{\prime}+\left(4f^{\prime}(r+2\xi^{(1)}\phi^{\prime})+f(4-\zeta r^{2}\phi^{\prime 2})\right)^{2}}\right).$

(14)

For the domain of the radial coordinate, we consider $r>r_{h}$ with $r_{h}$ ($>0$) being the position of the event (outermost) horizon of the spacetime. We assume that $f>0$, $f^{\prime}>0$, $h>0$, $h^{\prime}>0$ for $r>r_{h}$, and in the limit of $r\to r_{h}$ $f\to 0$, $h\to 0$, $\frac{f}{h}\to{\rm constant}$, and $\phi$ and $\phi^{\prime}$ are regular. Then, in the limit of $r\to r_{h}$, Eq. (12) becomes

$\displaystyle\left.\frac{f}{h}\right|_{r\to r_{h}}=\left.\left(r+2\xi^{(1)}\phi^{\prime}\right)f^{\prime}\right|_{r\to r_{h}},$

(15)

which can be obtained only from the $(+)$-branch of Eq. (14), provided that $r+2\xi^{(1)}\phi^{\prime}|_{r\to r_{h}}>0$. Thus, in the rest of the paper we only focus on the $(+)$-branch. Obviously, Eq. (15) implies that

$\displaystyle\left.\left(r+2\xi^{(1)}\phi^{\prime}\right)h^{\prime}\right|_{r\to r_{h}}=1.$

(16)

Assuming $r+2\xi^{(1)}\phi^{\prime}|_{r\to r_{h}}>0$ and substituting the $(+)-$branch of Eq. (14) into Eqs. (11) and (13), we obtain a set of the equations which are quasi-linear for $f^{\prime\prime}(r)$ and $\phi^{\prime\prime}(r)$, and, by rearranging them, we obtain the evolution equations for $f(r)$ and $\phi(r)$ with respect to $r$:

$\displaystyle f^{\prime\prime}(r)=F_{f}\left[f,f^{\prime},\phi,\phi^{\prime};r\right],\qquad\phi^{\prime\prime}(r)=F_{\phi}\left[f,f^{\prime},\phi,\phi^{\prime};r\right],$

(17)

where $F_{f}$ and $F_{\phi}$ are the nonlinear combinations of the given variables. After integrating Eq. (17) for $f(r)$ and $\phi(r)$ with the given boundary conditions near the event horizon and substituting them into Eq. (14), $h(r)$ can be obtained numerically.

We assume that $\xi(\phi)$ is a $Z_{2}$-symmetric function across $\phi=0$, $\xi(-\phi)=\xi(\phi)$, and thus satisfies $\xi^{(1)}(0)=\xi^{(3)}(0)=\xi^{(5)}(0)=\cdots=\xi^{(2i+1)}(0)=\cdots=0$ (where $i=3,5,7,\cdots$). Since the GB invariant is topological invariant, we may also set $\xi(0)=0$ without loss of generality, and hence the most general $Z_{2}$-symmetric coupling function which satisfies our requirement is given by

$\displaystyle\xi(\phi)=\frac{\eta}{8}\phi^{2}+\sum_{i=2}^{\infty}\alpha_{2i}\phi^{2i},$

(18)

where $\alpha_{4,6,8,\cdots}$ are the constant coefficients for the higher-order terms. In order to realize the tachyonic instability of the Schwarzschild solution against the radial perturbations, we require $\eta>0$.

Near the horizon $r=r_{h}$, the solution can be expanded as Silva et al. (2018); Doneva and Yazadjiev (2018a); Antoniou et al. (2018); Minamitsuji and Ikeda (2019); Silva et al. (2019)

	$\displaystyle f(r)$	$\displaystyle=$	$\displaystyle f_{1}\left(r-r_{h}\right)+{\cal O}\left((r-r_{h})^{2}\right),$		(19)
	$\displaystyle h(r)$	$\displaystyle=$	$\displaystyle\frac{r_{h}\zeta}{48\xi^{(1)}(\phi_{0})^{2}}\left(r_{h}^{2}-\sqrt{r_{h}^{4}-\frac{96\xi^{(1)}(\phi_{0})^{2}}{\zeta}}\right)\left(r-r_{h}\right)+{\cal O}\left((r-r_{h})^{2}\right),$		(20)
	$\displaystyle\phi(r)$	$\displaystyle=$	$\displaystyle\phi_{0}\left[1-\frac{1}{4r_{h}\phi_{0}\xi^{(1)}(\phi_{0})}\left(r_{h}^{2}-\sqrt{r_{h}^{4}-\frac{96\xi^{(1)}(\phi_{0})^{2}}{\zeta}}\right)\left(r-r_{h}\right)\right]+{\cal O}\left((r-r_{h})^{2}\right),$		(21)

where $\phi_{0}$ is the value of the scalar field at the event horizon, $f_{1}$ is an arbitrary constant representing the time-reparametrization invariance in the static and spherically symmetric spacetimes, and we omit to show the ${\cal O}\left((r-r_{h})^{2}\right)$ terms explicitly. As a consistency check, one can easily confirm that the above assumption $r+2\xi^{(1)}\phi^{\prime}|_{r\to r_{h}}>0$ and the relation (16) are satisfied, using the above expansions. On the other hand, in the large-distance limit $r\gg r_{h}$, we obtain the asymptotic solutions

	$\displaystyle\frac{f}{f_{\infty}}$	$\displaystyle=1-\frac{2M}{r}+\frac{\zeta MQ^{2}}{12r^{3}}+\frac{MQ}{6r^{4}}\left(\zeta MQ+24\xi^{(1)}(\phi_{\infty})\right)+{\cal O}\left(\frac{1}{r^{5}}\right),$
	$\displaystyle h$	$\displaystyle=1-\frac{2M}{r}+\frac{\zeta Q^{2}}{4r^{2}}+\frac{\zeta MQ^{2}}{4r^{3}}+\frac{MQ}{3r^{4}}\left(\zeta MQ+24\xi^{(1)}(\phi_{\infty})\right)+{\cal O}\left(\frac{1}{r^{5}}\right),$
	$\displaystyle\phi$	$\displaystyle=\phi_{\infty}+\frac{Q}{r}+\frac{MQ}{r^{2}}+\frac{1}{r^{3}}\left(\frac{4M^{2}Q}{3}-\frac{\zeta Q^{3}}{24}\right)+\frac{M}{6r^{4}}\left(12M^{2}Q-\zeta Q^{3}-\frac{24M\xi^{(1)}(\phi_{\infty})}{\zeta}\right)+{\cal O}\left(\frac{1}{r^{5}}\right),$		(22)

where $\phi_{\infty}:=\phi(r\to\infty)$ is the asymptotic value of the scalar field, and $M$ and $Q$ are the mass and the scalar charge, respectively. We note that the constant $f_{\infty}$ ($>0$) also represents the time-reparametrization invariance in the static and spherically symmetric spacetimes, and may be set to unity after the proper rescaling of the time coordinate. From Eq. (II.2), $M$ and $Q$ can be numerically evaluated as

$\displaystyle M=\frac{r}{2}\left(1-h\right)\Big{|}_{r\to\infty},\qquad Q=-r^{2}\phi^{\prime}(r)\Big{|}_{r\to\infty}.$

(23)

The scalarized solutions connect two values of the scalar field: $\phi_{0}\neq 0$ near the horizon $r=r_{h}$ and

$\displaystyle\phi_{\infty}=0,$

(24)

in the large-distance limit.

We focus on the nodeless scalarized solution where the scalar field $\phi$ monotonically approaches 0 from a nonzero value on the horizon, which is known as the fundamental solution and only the stable solution against the radial perturbations Blazquez-Salcedo et al. (2018); Minamitsuji and Ikeda (2019). In order to satisfy the boundary condition (24), we require that $\phi$ monotonically approaches $0$ from $\phi_{0}\neq 0$ on the horizon. From Eq. (21), we then have to impose $\phi_{0}\xi^{(1)}(\phi_{0})>0$ near the horizon. For the $Z_{2}$-symmetric coupling (18), without loss of generality we may assume that

$\displaystyle\phi_{0}>0,\qquad\xi^{(1)}(\phi_{0})>0.$

(25)

Under the boundary condition (24) for the scalarized solutions, with the use of the properties $\xi^{(1)}(0)=\xi^{(3)}(0)=\xi^{(5)}(0)=\cdots=\xi^{(2i+1)}(0)=\cdots=0$ (where $i=3,4,5,\cdots$), Eq. (II.2) reduces to

	$\displaystyle\frac{f}{f_{\infty}}$	$\displaystyle=1-\frac{2M}{r}+\frac{\zeta MQ^{2}}{12r^{3}}+\frac{\zeta M^{2}Q^{2}}{6r^{4}}+{\cal O}\left(\frac{1}{r^{5}}\right),$
	$\displaystyle h$	$\displaystyle=1-\frac{2M}{r}+\frac{\zeta Q^{2}}{4r^{2}}+\frac{\zeta MQ^{2}}{4r^{3}}+\frac{\zeta M^{2}Q^{2}}{3r^{4}}+{\cal O}\left(\frac{1}{r^{5}}\right),$
	$\displaystyle\phi$	$\displaystyle=\frac{Q}{r}+\frac{MQ}{r^{2}}+\frac{1}{r^{3}}\left(\frac{4M^{2}Q}{3}-\frac{\zeta Q^{3}}{24}\right)+\frac{M}{6r^{4}}\left(12M^{2}Q-\zeta Q^{3}\right)+{\cal O}\left(\frac{1}{r^{5}}\right).$		(26)

As we will see later, in order to check whether one of the stability conditions (36) introduced later is satisfied or not, we need terms higher-order in the power of $r^{-1}$, which are omitted in Eq. (II.2). For the stability analysis of the scalarized solution (see Sec. IV) we need to expand the background solution at least up to ${\cal O}(r^{-8})$, and for the stability analysis of the nonlinearly scalarized solutions (see Sec. V) we need to do that at least up to ${\cal O}(r^{-12})$. However, since the expression of these higher-order corrections to Eq. (II.2) are quite involved, we do not show them explicitly here.

In the presence of matter fields, scalarization of NSs could also be realized for the same type of the coupling functions Doneva and Yazadjiev (2018b). In contrast to the BH scalarization in the static and spherically symmetric spacetime, since the GB could change the sign in the presence of matter fields, the NS scalarization could occur for both signs of the coupling terms. Assuming the same asymptotic value of the scalar field as Eq. (24), the metric and the scalar field in the exterior solutions for scalarized NSs can also be described by Eq. (II.2). In this case, the values of the mass and scalar charge, $M$ and $Q$, are determined via the matching with the interior solutions with matter fields at the surface of a NS and the use of Eq. (23) in the asymptotic limit.

The linear stability against the odd-parity perturbations is ensured under the following three conditions Kobayashi et al. (2012):

	$\displaystyle{\cal F}$	$\displaystyle:=$	$\displaystyle 2G_{4}+h\phi^{\prime 2}G_{5,\phi}-h\phi^{\prime 2}\left(\frac{1}{2}h^{\prime}\phi^{\prime}+h\phi^{\prime\prime}\right)G_{5,X}>0\,,$		(27)
	$\displaystyle{\cal G}$	$\displaystyle:=$	$\displaystyle 2G_{4}+2h\phi^{\prime 2}G_{4,X}-h\phi^{\prime 2}\left(G_{5,\phi}+{\frac{f^{\prime}h\phi^{\prime}G_{5,X}}{2f}}\right)>0\,,$		(28)
	$\displaystyle{\cal H}$	$\displaystyle:=$	$\displaystyle 2G_{4}+2h\phi^{\prime 2}G_{4,X}-h\phi^{\prime 2}G_{5,\phi}-\frac{h^{2}\phi^{\prime 3}G_{5,X}}{r}>0\,.$		(29)

The squared propagation speeds of the odd-parity perturbations along the radial and angular directions are given, respectively, by

$\displaystyle c_{r,{\rm odd}}^{2}=\frac{{\cal G}}{{\cal F}}\,,\qquad c_{\Omega,{\rm odd}}^{2}=\frac{{\cal G}}{{\cal H}}\,.$

(30)

Thus, if all the conditions (27)-(29) are satisfied, all the sound speeds in Eq. (30) are positive.

In the even-parity perturbations, the kinetic term of the tensorial polarization has the correct sign for (28), and then that for the scalar field has the correct sign Kobayashi et al. (2014), if the following condition is satisfied

$\displaystyle{\cal K}:=2{\cal P}_{1}-{\cal F}>0\,,$

(31)

with

$\displaystyle{\cal P}_{1}:=\frac{h\mu}{2fr^{2}{\cal H}^{2}}\left(\frac{fr^{4}{\cal H}^{4}}{\mu^{2}h}\right)^{\prime}\,,\qquad\mu:=\frac{2(\phi^{\prime}a_{1}+r\sqrt{fh}{\cal H})}{\sqrt{fh}}\,,$

(32)

where $a_{1}$ is given in Appendix A. In the limit of high frequencies, the conditions for the absence of Laplacian instabilities of the even-parity tensorial polarization $\psi$ and the scalar field polarization $\delta\phi$ along the radial direction are given, respectively, by

	$\displaystyle c_{r1,{\rm even}}^{2}$	$\displaystyle=$	$\displaystyle\frac{\mathcal{G}}{\mathcal{F}}>0\,,$		(33)
	$\displaystyle c_{r2,{\rm even}}^{2}$	$\displaystyle=$	$\displaystyle\frac{2\phi^{\prime}[4r^{2}(fh)^{3/2}{\cal H}c_{4}(2\phi^{\prime}a_{1}+r\sqrt{fh}\,{\cal H})-2a_{1}^{2}f^{3/2}\sqrt{h}\phi^{\prime}{\cal G}+(a_{1}f^{\prime}+2c_{2}f)r^{2}fh{\cal H}^{2}]}{f^{5/2}h^{3/2}(2{\cal P}_{1}-{\cal F})\mu^{2}}>0\,,$		(34)

where $c_{2}$ and $c_{4}$ are presented in Appendix A. Since $c_{r1,{\rm even}}^{2}$ is the same as $c_{r,{\rm odd}}^{2}$, only the second propagation speed squared $c_{r2,{\rm even}}^{2}$ provides an additional stability condition. We note that for the monopole mode $\ell=0$ there is no propagation for the gravitational perturbation, while the scalar-field perturbation $\delta\phi$ propagates with the same radial velocity as Eq. (34). We also note that for the dipole mode $\ell=1$ there is only one gauge DOF for fixing $\delta\phi=0$, under which the gravitational perturbation propagates with the same radial speed squared as Eq. (34).

We then turn to the linear stability conditions against the propagation in the angular directions. In the limit of large multipoles $\ell\gg 1$, the conditions associated with the squared angular propagation speeds in the even-parity perturbations are Kase and Tsujikawa (2022); Minamitsuji et al. (2022a, b)

$\displaystyle c_{\Omega\pm}^{2}=-B_{1}\pm\sqrt{B_{1}^{2}-B_{2}}>0\,,$

(35)

where we present the explicit form of $B_{1}$ and $B_{2}$ in Appendix A. These conditions are satisfied if and only if

$\displaystyle B_{1}^{2}\geq B_{2}>0\quad{\rm and}\quad B_{1}<0\,.$

(36)

First, we consider the quartic-order coupling function discussed in Ref. Minamitsuji and Ikeda (2019); Silva et al. (2019)

$\displaystyle\xi(\phi)=\frac{\eta}{8}\left(\phi^{2}+\alpha\phi^{4}\right),$

(37)

where $\eta$ and $\alpha$ are constants. We require that $\eta>0$ and $\alpha<0$, so that the Schwarzschild solution suffers from the tachyonic instability and the scalarized BHs are stable against the radial perturbations Minamitsuji and Ikeda (2019); Silva et al. (2019), respectively.

In the limit of $r\to\infty$, under the boundary condition (24), the functions $\mathcal{F}$, $\mathcal{G}$, and $\mathcal{H}$ defined in Eqs. (27)-(29) can be expanded as

$\displaystyle\mathcal{F}$

$\displaystyle=$

$\displaystyle\frac{1}{\kappa^{2}}\left(1-\frac{3Q^{2}\eta}{r^{4}}\right)+{\cal O}\left(\frac{1}{r^{6}}\right),\quad\mathcal{G}=\frac{1}{\kappa^{2}}\left(1+\frac{MQ^{2}\eta}{r^{5}}\right)+{\cal O}\left(\frac{1}{r^{6}}\right),\quad\mathcal{H}=\frac{1}{\kappa^{2}}\left(1+\frac{Q^{2}\eta}{r^{4}}\right)+{\cal O}\left(\frac{1}{r^{6}}\right),$

(38)

which are always positive at the leading order. The sound speeds for the radial and angular propagations in the odd-parity perturbations (30) and that for the radial propagation of the metric perturbations in the even-parity perturbations (33) coincide with the speed of light at the leading order, with the corrections of ${\cal O}(r^{-4})$. The function $\mathcal{K}$ defined by Eq. (31) can be expanded as

$\displaystyle\mathcal{K}=\frac{\zeta Q^{2}}{\kappa^{2}r^{2}}\left(\frac{1}{4}+\frac{M}{r}\right)+{\cal O}\left(\frac{1}{r^{4}}\right),$

(39)

which is also positive for the correct sign of the kinetic term $\zeta>0$. The sound speed of the radial perturbations of the scalar field in the even-parity perturbations $c_{r2,{\rm even}}^{2}$, which can be evaluated via Eq. (34), also coincides with the speed of light at the leading order, with the corrections of ${\cal O}(r^{-9})$. Thus, in the large $r$ region the scalarized BHs are linearly stable against all types of propagations in the odd-parity perturbations and against the radial propagations in the even-parity perturbations.

We now check the angular sound speeds of the even-parity perturbations. The functions $B_{1}$ and $B_{2}$ can be expanded as

	$\displaystyle B_{1}$	$\displaystyle=$	$\displaystyle-1+\frac{Q^{2}\eta}{2r^{4}}+\frac{Q^{2}\left(-20M^{2}+Q^{2}(24\alpha+\zeta)\right)\eta}{24r^{6}}+\frac{MQ^{2}\eta\left(-16M^{2}\zeta+Q^{2}\zeta(12\alpha+\zeta)+12\eta\right)}{6\zeta r^{7}}+{\cal O}\left(\frac{1}{r^{8}}\right),$		(40)
	$\displaystyle B_{2}$	$\displaystyle=$	$\displaystyle 1-\frac{Q^{2}\eta}{r^{4}}-\frac{Q^{2}\left(-20M^{2}+Q^{2}(24\alpha+\zeta)\right)\eta}{12r^{6}}+\frac{MQ^{2}\eta\left(16M^{2}\zeta-Q^{2}\zeta(12\alpha+\zeta)+24\eta\right)}{3\zeta r^{7}}+{\cal O}\left(\frac{1}{r^{8}}\right),$		(41)

leading to

$\displaystyle B_{1}^{2}-B_{2}=-\frac{12MQ^{2}\eta^{2}}{\zeta r^{7}}+{\cal O}\left(\frac{1}{r^{8}}\right),$

(42)

which is negative at the leading order and the first condition of Eq. (36) is not satisfied, since $M>0$ and $\zeta>0$ for the correct sign of the scalar kinetic term. We also note that the above result (42) is independent of the sign of $\eta$ and $Q$. In addition, since the leading term in Eq. (42) does not depend on the coefficient of the quartic-order term $\alpha$, we expect that the same leading behavior should be obtained for other nonlinear coupling functions.

To confirm our expectation in the previous subsection, we consider the exponential coupling function discussed originally in Ref. Doneva and Yazadjiev (2018a)

$\displaystyle\xi(\phi)=\frac{\eta}{8\beta}\left(1-e^{-\beta\phi^{2}}\right),$

(43)

which has also been employed in the literature Doneva and Yazadjiev (2018a); Blazquez-Salcedo et al. (2018); Cunha et al. (2019); Doneva and Yazadjiev (2021); East and Ripley (2021); Doneva et al. (2022b). Again, in order to obtain BH solutions which are linearly stable against the radial perturbations, we require that $\eta>0$ and $\beta>0$. In the limit of the spatial infinity, $r\to\infty$, under the boundary condition (24), the functions of $\mathcal{F}$, $\mathcal{G}$, and $\mathcal{H}$ can be expanded as Eq. (38), which are always positive. The function $\mathcal{K}$ can be expanded as Eq. (39), which is also positive at the leading order. The sound speeds for the radial and angular propagations in the odd-parity perturbations (30) and that for the radial propagation of the metric perturbations in the even-parity perturbations (33) coincide with the speed of light at the leading order, with the corrections of ${\cal O}(r^{-4})$. The sound speed of the radial propagation of the scalar field in the even-parity perturbations $c_{r2,{\rm even}}^{2}$ can be evaluated via Eq. (34) and also coincides with the speed of light at the leading order, with the corrections of ${\cal O}(r^{-9})$. Thus, in the large $r$ region the scalarized BHs are linearly stable against the odd-parity perturbations and against the radial propagations in the even-parity perturbations.

The functions $B_{1}$ and $B_{2}$ can be expanded as

	$\displaystyle B_{1}$	$\displaystyle=$	$\displaystyle-1+\frac{Q^{2}\eta}{2r^{4}}+\frac{Q^{2}\left(-20M^{2}+Q^{2}(-12\beta+\zeta)\right)\eta}{24r^{6}}+\frac{MQ^{2}\eta\left(-16M^{2}\zeta+Q^{2}\zeta(-6\beta+\zeta)+12\eta\right)}{6\zeta r^{7}}+{\cal O}\left(\frac{1}{r^{8}}\right),$		(44)
	$\displaystyle B_{2}$	$\displaystyle=$	$\displaystyle 1-\frac{Q^{2}\eta}{r^{4}}+\frac{Q^{2}\left(20M^{2}+Q^{2}(12\beta-\zeta)\right)\eta}{12r^{6}}+\frac{MQ^{2}\eta\left(16M^{2}\zeta+Q^{2}\zeta(6\beta-\zeta)+24\eta\right)}{3\zeta r^{7}}+{\cal O}\left(\frac{1}{r^{8}}\right),$		(45)

which also lead to the same leading behavior of $B_{1}^{2}-B_{2}$ as Eq. (42), and again we find that the first condition of Eq. (36) is not satisfied.

Ref. Doneva and Yazadjiev (2018b) showed that the same coupling function as Eq. (43) could realize scalarization of static and spherically symmetric NSs. In contrast to the case of BH scalarization, since in the static and spherically symmetric spacetimes the GB invariant could change the sign in the presence of matter fields, NS scalarization could occur for both signs of the parameter $\eta$. Since the metric and the scalar field in the exterior solutions for scalarized NSs are also described by Eq. (II.2), our results in the section can also be applied to the case of NS scalarization. As the leading-order term in Eq. (42) is irrespective of the sign of $\eta$, the first condition of Eq. (36) is violated also for NS scalarization with any matter equation of state. Thus, our analysis in this section should exclude both BH and NS scalarization models induced by the GB coupling function of the form (43).

Along the same analysis for more general couplings (18), we obtain the same leading behavior as Eq. (42). Thus, the instability does not depend on the higher-order structure of the coupling function $\xi(\phi)$. The result that the gradient instability in the even-parity perturbations appears in the angular directions, irrespective of the detailed structure of $\xi(\phi)$, implies that in the large-$\ell$ limit the onset of this instability would take place in the vicinity of the bifurcation point of the scalarized branch from the Schwarzschild branch in the mass-charge diagram, i.e., on the axis of $Q=0$, whose position is also irrespective of the higher-order terms in $\xi(\phi)$.

It would be very interesting if one can obtain universal constraints on the generic coupling functions $\xi(\phi)$ from theoretical considerations. For example, it is known that so-called positivity bounds Adams et al. (2006) put non-trivial constraints on low-energy effective field theories by assuming the existence of local, causal, unitary, and Lorentz-invariant ultraviolet (UV) completions. However, the positivity bounds in the presence of gravity are rather subtle and still at the stage of development, in particular require some additional assumptions about unknown behaviors of the UV completion of gravity (see e.g. Ref. Aoki et al. (2021a)). Moreover, even without inclusion of gravity the naive positivity bounds may be violated around Lorentz-violating backgrounds Aoki et al. (2021b). Nonetheless, some trials to consider positivity bounds in the context of the Einstein-scalar-GB theories have been made (see e.g. Ref. Herrero-Valea (2022)). It is certainly worthwhile developing our understanding of positivity bounds in the presence of gravity around Lorentz-violating backgrounds such as BHs.