# Even-parity stability of hairy black holes in $U(1)$ gauge-invariant
scalar-vector-tensor theories

Chao Zhang1,2 and Ryotaro Kase3 1Basic Research Center for Energy
Interdisciplinary, College of Science, China University of Petroleum-Beijing,
Beijing 102249, China.
2Beijing Key Laboratory of Optical Detection Technology for Oil and Gas, China
University of Petroleum-Beijing, Beijing, 102249, China.
3 Department of Physics, Faculty of Science, Tokyo University of Science, 1-3,
Kagurazaka, Shinjuku-ku, Tokyo 162-8601, Japan.

###### Abstract

The $U(1)$ gauge-invariant scalar-vector-tensor theories, which catches five
degrees of freedom, are valuable for its implications to inflation problems,
generation of primordial magnetic fields, new black hole (BH) and neutron star
solutions, etc. In this paper, we derive conditions for the absence of ghosts
and Laplacian instabilities of nontrivial BH solutions dressed with scalar
hair against both odd- and even-parity perturbations on top of the static and
spherically symmetric background in the most general $U(1)$ gauge-invariant
scalar-vector-tensor theories with second-order equations of motion. In
addition to some general discussions, several typical concrete models are
investigated. Specially, we show that the stability against even-parity
perturbations is ensured outside the event horizon under certain constraints
to these models. This is a crucial step to check the self-consistency of the
theories and to shed light on the physically accessible models of such
theories for future studies.

##  I Introduction

Over the past decades, general relativity (GR) has passed all the experimental
tests with flying colors (see, e.g., Ref. Özel and Freire (2016)). From the
theoretical point of view, it is quite an elegant, robust and systematic
theory to describe the gravity. Nonetheless, there are still many crucial open
questions left by the framework of GR. For instance, the accelerated expansion
of the universe, the origin and inherence of dark matter/energy, the
possibility to construct a quantum theory of gravity, etc., are some
outstanding ones Debono and Smoot (2016). To solve these problems and approach
to the nature of gravity, many modified gravitational theories were
constructed by introducing additional interactions and fields. For the
purposes of justifying and confining those modified gravitational theories,
they have to be subjected to experimental tests, including those from the
solar system Will (2014) as well as strong-field regimes, e.g., the
gravitational waves (GWs) led by compact celestial bodies Tattersall  _et al._
(2018); Barack  _et al._ (2019); Baibhav  _et al._ (2021).

The detection of the first GW from the coalescence of two massive black holes
(BHs) by advanced LIGO/Virgo marked the beginning of a new era — the GW
astronomy Abbott  _et al._ (2016). Following this observation, more than 90 GW
events have been identified by the LIGO/Virgo/KAGRA (LVK) scientific
collaborations (see, e.g., Refs. Abbott  _et al._ (2019, 2021a, 2020, 2021b)).
In the future, more ground- and space-based GW detectors will be constructed
Moore  _et al._ (2015), which will enable us to probe signals with a much
wider frequency band and larger distances. As a result, more types of GW
sources will be realized, including those from a final remnant BH Shi  _et
al._ (2019). BH is one of the most mysterious phenomena in the universe. The
existence of BHs provides us a perfect way to test gravitational effects under
extremely strong gravity as we pursued and mentioned above. On the other hand,
from the theoretical point of view, BHs are also unique laboratories to test
the deviation of modified theories of gravity from GR. The outbreak of
interest on BHs has further gained momenta after the detection of the shadow
of the supermassive compact objects in the center of galaxy M87 and
Sagittarius A* (Sgr A*) in the center of the Milky Way galaxy, which are the
most likely black hole candidates, with the Event Horizon Telescope (EHT)
Akiyama  _et al._ (2019a, b, 2021, 2022).

The development of the GW astronomy as well as the interest on BHs have
triggered the interest in the quasi-normal mode (QNM) of BHs, as GWs emitted
in the ringdown phase can be considered as the linear combination of these
individual modes Berti  _et al._ (2018a). From the classical point of view,
QNMs are eigenmodes of dissipative systems. The information contained in QNMs
provides the keys to revealing whether BHs are ubiquitous in our universe, and
more importantly whether GR is the correct theory to describe gravity even in
the strong field regime Berti  _et al._ (2018b). In addition to the
observational purposes, QNM is also an important indicator of the stability of
a specific spacetime Lin  _et al._ (2016); Zhang  _et al._ (2023a). Under
certain circumstance, the results from QNMs will exhibit a manifest
consistency with those from the Lagrangian-based stability analysis Tsujikawa
_et al._ (2021); Zhang  _et al._ (2023b).

In GR, according to the no-hair theorem, an isolated and stationary BH is
completely characterized by only three quantities, mass, angular momentum, and
electric charge. Astrophysically, we expect BHs to be electrically neutral, so
they are uniquely described by the Kerr solution. Nonetheless, in theories
that beyond GR, the existence of additional degrees of freedom (DOFs) can give
rise to new hairs to the field configuration and spacetime metric Heisenberg
_et al._ (2018). The theories containing a scalar field $\phi$ coupled to
gravity besides two tensor polarizations arising from the gravity sector are
dubbed scalar-tensor theories. In particular, Horndeski constructed most
general scalar-tensor theories with second-order equations of motion Horndeski
(1974); Deffayet  _et al._ (2011); Kobayashi  _et al._ (2011). On the other
hand, for a vector field coupled to gravity, it is known that generalized
Proca theories are the most general vector-tensor theories with second-order
equations of motion Heisenberg (2014); Tasinato (2014a, b) (see also Refs.
Errasti Díez  _et al._ (2020a, b)). These two important classes of field
theories, Horndeski and generalized Proca theories, can be unified in the
framework of scalar-vector-tensor (SVT) theories with second-order equations
of motion Heisenberg (2018). The SVT theories can be classified into two cases
depending on whether they respect the $U(1)$ gauge symmetry or not. In the
presence of $U(1)$ gauge symmetry the longitudinal component of a vector field
vanishes, so that the propagating DOFs are five in total (one scalar, two
transverse vectors, two tensor polarizations). The breaking of $U(1)$ gauge
symmetry leads to the propagation of the longitudinal scalar besides the five
DOFs. Although a canonical scalar field in the Einstein-Maxwell theories
cannot have non-trivial profile under the assumption of static and spherically
symmetric configuration, i.e., the no-hair theorem Chase (1970); Bekenstein
(1972, 1996) holds in the theories, this is not the case in SVT theories due
to the existence of the coupling between scalar and vector DOFs. Indeed, one
can construct a hairy, static and spherically symmetric BH solutions in $U(1)$
gauge-invariant SVT theories in which the scalar field can possess nontrivial
profileHeisenberg and Tsujikawa (2018); Ikeda  _et al._ (2019).

In this paper, we study the stability of static and spherically symmetric BHs
in the $U(1)$ gauge-invariant (GI) SVT theory. The one against the odd-parity
perturbations has been studied in Ref. Heisenberg  _et al._ (2018). This paper
is a successor that turns to the even-parity sector Thompson  _et al._ (2017);
Kase and Tsujikawa (2021, 2022); Liu  _et al._ (2023), which completes this
kind of Lagrangian-based stability analysis. By demanding the self-consistency
of the theory, Ref. Heisenberg  _et al._ (2018) has already brought some
constraints to the phase space of the theory-dependent coupling parameters for
certain models of the $U(1)$ GI SVT theory. As will be seen later, in
combining with the even-parity analysis, the corresponding phase space will be
further confined.

The rest of the paper is organized as following: Sec. II provides the
necessary fundamental information of the $U(1)$ GI SVT theory. The
corresponding background field equations are discussed there. After that, we
quickly review the stability analysis against the odd-parity perturbations in
Sec. III. Taking advantages of the results led by the odd-parity sector, we
further run the stability analysis for the even-parity case in Sec. IV. Notice
that, the analysis is divided into 3 parts according to $l\geq 2$, $l=0$ and
$l=1$. We apply our general stability conditions to the three typical models
in Sec. V for the background solutions studied in Ref. Heisenberg and
Tsujikawa (2018). Finally, some of the concluding remarks will be given in
Sec. VI.

As a usual treatment, in the following we shall set the speed of light as well
as the reduced Planck constant to one, viz., $c=\hbar=1$111Notice that, after
this unit selection, there is still one DOF left for the unit system of
$[L,M,T]$. As an example, one can further set the radii of metric horizon to
the Planck mass ($r_{h}=M_{\text{pl}}$) to fix the unit system.. All the Greek
letters in indices run from 0 to 3. Other usages of indices will be explained
when it is necessary. The whole paper is working under the signature $(-+++)$.

##  II Background equations in $U(1)$ gauge-invariant SVT theories

We consider the $U(1)$ gauge-invariant (GI) scalar-vector-tensor (SVT)
theories described by the action

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

where $g$ is a determinant of the metric tensor $g_{\mu\nu}$, $M_{\rm
pl}^{2}R/2$ is the Einstein-Hilbert term composed of the reduced Planck mass
$M_{\rm pl}$ associated with Newton’s gravitational constant $G$ as $M_{\rm
pl}=1/\sqrt{8\pi G}$ and $R$ denotes the Ricci scalar. The Lagrangians
$\mathcal{L}^{i}_{\rm SVT}$ with $i=2,3,4$ representing the $U(1)$ GI SVT
interactions Heisenberg (2018) between a scalar field $\phi$ and a vector
field $A_{\mu}$ are given by

$\displaystyle\mathcal{L}^{2}_{\rm SVT}$ $\displaystyle=$ $\displaystyle
f_{2}(\phi,X,F,\tilde{F},Y)\,,$ (2) $\displaystyle\mathcal{L}^{3}_{\rm SVT}$
$\displaystyle=$
$\displaystyle\left[f_{3}(\phi,X)g_{\rho\sigma}+\tilde{f}_{3}(\phi,X)\nabla_{\rho}\phi\nabla_{\sigma}\phi\right]\tilde{F}^{\mu\rho}\tilde{F}^{\nu\sigma}\nabla_{\mu}\nabla_{\nu}\phi\,,$
(3) $\displaystyle\mathcal{L}^{4}_{\rm SVT}$ $\displaystyle=$ $\displaystyle
f_{4}(\phi,X)L^{\mu\nu\alpha\beta}F_{\mu\nu}F_{\alpha\beta}+\left[\frac{1}{2}f_{4,X}(\phi,X)+\tilde{f}_{4}(\phi)\right]\tilde{F}^{\mu\nu}\tilde{F}^{\alpha\beta}\nabla_{\mu}\nabla_{\alpha}\phi\nabla_{\nu}\nabla_{\beta}\phi\,,$
(4)

where $\nabla_{\mu}$ is the covariant derivative operator. The function
$f_{2}$ depends on $\phi$ and the following quantities,

$X\equiv-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi\,,\qquad
F\equiv-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}\,,\qquad\tilde{F}\equiv-\frac{1}{4}F_{\mu\nu}\tilde{F}^{\mu\nu}\,,\qquad
Y\equiv\nabla_{\mu}\phi\nabla_{\nu}\phi F^{\mu\alpha}{F^{\nu}}_{\alpha}\,,$
(5)

where

$F_{\mu\nu}\equiv\nabla_{\mu}A_{\nu}-\nabla_{\nu}A_{\mu}\,,\qquad\tilde{F}^{\mu\nu}\equiv\frac{1}{2}\mathcal{E}^{\mu\nu\alpha\beta}F_{\alpha\beta}\,.$
(6)

Here, the anti-symmetric Levi-Civita tensor $\mathcal{E}^{\mu\nu\alpha\beta}$
satisfies the normalization
$\mathcal{E}^{\mu\nu\alpha\beta}\mathcal{E}_{\mu\nu\alpha\beta}=-4!$. We
denote the derivative of $f_{2}$ as $f_{2,Z}=\partial f_{2}/\partial Z$ where
$Z$ is any of $\phi$, $X$, $F$, $\tilde{F}$, and $Y$. Meanwhile,
$f_{3},\tilde{f}_{3},f_{4}$ are functions of $\phi,X$ with the same notation
such as $f_{3,\phi}=\partial f_{3}/\partial\phi$, $f_{3,X}=\partial
f_{3}/\partial X$, etc., and $\tilde{f}_{4}$ depends on $\phi$ alone as seen
from Eq. (4). The double dual Riemann tensor $L^{\mu\nu\alpha\beta}$ is
defined by

$L^{\mu\nu\alpha\beta}\equiv\frac{1}{4}\mathcal{E}^{\mu\nu\rho\sigma}\mathcal{E}^{\alpha\beta\gamma\delta}R_{\rho\sigma\gamma\delta}\,,$
(7)

where $R_{\rho\sigma\gamma\delta}$ is the Riemann tensor.

In this paper, we study the stability of the hairy BH solutions in the $U(1)$
GI theories studied in Ref. Heisenberg and Tsujikawa (2018) on top of the
static and spherically symmetric background given by the line element [under
the Boyer-Lindquist coordinate $(t,r,\theta,\varphi)$]

${\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)\,,$ (8)

where $f$ and $h$ depend on the radial coordinate $r$. According to the
underlying symmetry of the spacetime, we consider the scalar field $\phi$
depending on $r$ alone at the level of background, such that

$\phi=\phi(r)\,.$ (9)

Similarly, the background components of $A_{\mu}$ are given as

$A_{\mu}=(A_{0}(r),0,0,0)\,,$ (10)

where the radial component is absent due to the $U(1)$ gauge-invariance
Heisenberg  _et al._ (2018). We note that one can introduce the magnetic
charge $P$ by setting $A_{\varphi}=-P\cos\theta$ as in Refs. Lee and Weinberg
(1991); Fernandes  _et al._ (2019); Taniguchi  _et al._ (2024). However, we do
not include such term in this paper since we focus on the stability analysis
of the solutions given in Ref. Heisenberg and Tsujikawa (2018) in which the
magnetic charge is absent. Denoting the quantities $X$, $F$, $\tilde{F}$, and
$Y$ evaluated on the background with the overbar, they reduce to

$\overline{X}=-\frac{h}{2}\phi^{\prime
2}\,,\qquad\overline{F}=\frac{hA_{0}^{\prime
2}}{2f}\,,\qquad\overline{\tilde{F}}=0\,,\qquad\overline{Y}=4\overline{X}\,\overline{F}\,.$
(11)

Here, a prime in the superscript denotes the derivative with respect to $r$.
Since the dependence on $\tilde{F}$ and $Y$ in $f_{2}$ under a static and
spherically symmetric background either vanishes or can be expressed in terms
of $X$ and $F$, it can be omitted at the background level Heisenberg  _et al._
(2017); Kase  _et al._ (2018a). However, since the above relations hold only
at the background level, the dependence on $\tilde{F}$ and $Y$ in $f_{2}$ may
give rise to specific effect on the dynamics of the odd- and even-parity
perturbations. Thus, we keep the full dependence in $f_{2}$, i.e.,
$f_{2}=f_{2}(\phi,X,F,\tilde{F},Y)$, in this paper222In fact, as seen from
Sec. III, the quantity ${\tilde{F}}$ will hold the linear-order contributions
in the odd-party sector. In contrast, ${\tilde{F}}$ is a quantity with a
magnitude at most the third order of the gravitational perturbation in the
even-parity sector and hardly affects the linear perturbation calculations as
we will see in Sec. IV.. We omit the overbar in the following discussion
before stimulating any confusions.

By the variation of the action (1) with respect to $f,h,\phi,A_{0}$, we obtain
background equations of motion, respectively, as333We notice that the explicit
dependence of $f_{2}$ on $Y$ was not considered in the counterparts of Eqs.
(12)-(18) in Ref. Heisenberg and Tsujikawa (2018) by assuming that $Y$ has
been spelled out as $\overline{Y}=4\overline{X}\,\overline{F}$ at the
background as we can see in Eq. (11). In this paper, we explicitly include
such a dependence in $f_{2}$ to be consistent with the fact that the relation
$\overline{Y}=4\overline{X}\,\overline{F}$ does not hold at the perturbation
level. Thus, Eqs. (12)-(18) look slightly different in comparing to the
equations in Ref. Heisenberg and Tsujikawa (2018).

$\displaystyle{\cal E}_{00}$ $\displaystyle\equiv$ $\displaystyle M_{\rm
pl}^{2}rfh^{\prime}-\Big{[}M_{\rm
pl}^{2}f(1-h)+r^{2}\left\\{ff_{2}-hA_{0}^{\prime 2}(f_{2,F}-2h\phi^{\prime
2}f_{2,Y})\right\\}-2rh^{2}\phi^{\prime}A_{0}^{\prime 2}f_{3}$ (12)
$\displaystyle+hA_{0}^{\prime 2}\left\\{4(h-1)f_{4}-h^{2}\phi^{\prime
2}(f_{4,X}+2\tilde{f}_{4})\right\\}\Big{]}=0\,,$ $\displaystyle{\cal E}_{11}$
$\displaystyle\equiv$ $\displaystyle M_{\rm pl}^{2}rhf^{\prime}-\Big{[}M_{\rm
pl}^{2}f(1-h)+r^{2}\left\\{ff_{2}+fh\phi^{\prime 2}f_{2,X}-hA_{0}^{\prime
2}(f_{2,F}-4h\phi^{\prime
2}f_{2,Y})\right\\}-2rh^{2}\phi^{\prime}A_{0}^{\prime
2}\left(3f_{3}-h\phi^{\prime 2}f_{3,X}\right)$ (13)
$\displaystyle+hA_{0}^{\prime 2}\left\\{4(3h-1)f_{4}-h(9h-4)\phi^{\prime
2}f_{4,X}+h^{3}\phi^{\prime 4}f_{4,XX}-10h^{2}\phi^{\prime
2}\tilde{f}_{4}\right\\}\Big{]}=0\,,$ $\displaystyle{\cal E}_{\phi}$
$\displaystyle\equiv$ $\displaystyle J_{\phi}^{\prime}-{\cal P}_{\phi}=0\,,$
(14) $\displaystyle{\cal E}_{A}$ $\displaystyle\equiv$ $\displaystyle
J_{A}^{\prime}=0\,,$ (15)

where

$\displaystyle J_{\phi}$ $\displaystyle\equiv$
$\displaystyle-\sqrt{\frac{h}{f}}\left[r^{2}\left(ff_{2,X}+2hA_{0}^{\prime
2}f_{2,Y}\right)\phi^{\prime}-2hA_{0}^{\prime
2}\left(2h\tilde{f}_{4}+3hf_{4,X}-2f_{4,X}\right)\phi^{\prime}+2rh^{2}A_{0}^{\prime
2}f_{3,X}\phi^{\prime 2}\right.$ (16) $\displaystyle\left.+h^{3}A_{0}^{\prime
2}f_{4,XX}\phi^{\prime 3}-2rhA_{0}^{\prime 2}f_{3}\right]\,,$
$\displaystyle{\cal P}_{\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{fh}}\left[r^{2}ff_{2,\phi}+hA_{0}^{\prime
2}\left\\{4f_{4,\phi}+2h(r\phi^{\prime}f_{3,\phi}-2f_{4,\phi})+h^{2}\left(f_{4,X\phi}+2\tilde{f}_{4,\phi}\right)\phi^{\prime
2}\right\\}\right]\,,$ (17) $\displaystyle J_{A}$ $\displaystyle\equiv$
$\displaystyle\sqrt{\frac{h}{f}}A_{0}^{\prime}\left[r^{2}\left(f_{2,F}-2h\phi^{\prime
2}f_{2,Y}\right)+4rh\phi^{\prime}f_{3}+8(1-h)f_{4}+2h^{2}\phi^{\prime
2}\left(f_{4,X}+2\tilde{f}_{4}\right)\right]\,.$ (18)

From Eq. (15), the current $J_{A}$ of vector field is conserved by virtue of
the $U(1)$ gauge symmetry. In addition, if we demand the shift symmetry of the
scalar field, the current $J_{\phi}$ will also be conserved. In such a case,
this current vanishes at the horizon and hence it should be zero everywhere as
we will see in Sec. II. We note that the coupling $\tilde{f}_{3}$ never
appears in Eqs. (12)-(15) due to the underlying symmetry of background
spacetime.

##  III stability conditions against odd-parity perturbations

In this section, we revisit the stability conditions against odd-parity
perturbations in the $U(1)$ GI SVT theories which is studied in Ref.
Heisenberg  _et al._ (2018). We consider the perturbed metric
$g_{\mu\nu}=\bar{g}_{\mu\nu}+h_{\mu\nu}$ where $\bar{g}_{\mu\nu}$ is the
background metric defined by Eq. (8) and $h_{\mu\nu}$ is the perturbation
satisfying $|h_{\mu\nu}|\ll|\bar{g}_{\mu\nu}|$. The perturbations $h_{\mu\nu}$
can be expanded in terms of the spherical harmonics $Y_{lm}(\theta,\varphi)$
due to the underlying symmetry of the background spacetime. In doing so, the
perturbations are classified into odd- and even-parity modes where the former
changes the sign as $(-1)^{l+1}$ under the parity transformation
$(\theta,\varphi)\to(\pi-\theta,\varphi+\pi)$ and the latter as $(-1)^{l}$
Zerilli (1970). The components of $h_{\mu\nu}$ in the odd-parity perturbations
are expressed as follows

$\displaystyle h_{tt}=h_{tr}=h_{rr}=0\,,$ $\displaystyle
h_{ta}=\sum_{l,m}Q_{lm}(t,r)E_{ab}\nabla^{b}Y_{lm}(\theta,\varphi)\,,\qquad
h_{ra}=\sum_{l,m}W_{lm}(t,r)E_{ab}\nabla^{b}Y_{lm}(\theta,\varphi)\,,$
$\displaystyle
h_{ab}=\frac{1}{2}\sum_{l,m}U_{lm}(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]\,,$
(19)

where $a,b$ represent either $\theta$ or $\varphi$, $Q_{lm}$, $W_{lm}$, and
$U_{lm}$ are functions of $t$ and $r$. The tensor $E_{ab}$ is associated with
the anti-symmetric symbol $\varepsilon_{ab}$ satisfying
$\varepsilon_{\theta\varphi}=1$ as
$E_{ab}\equiv\sqrt{\gamma}\,\varepsilon_{ab}$, where
$\gamma\equiv\sin^{2}\theta$ is the determinant of the metric $\gamma_{ab}$
defined on the surface of two-dimension unit sphere. The scalar filed $\phi$
does not contribute to the odd-parity perturbations, i.e., it has
contributions only at the background level as long as we consider odd-parity
perturbations. The odd-parity vector field perturbations $\delta A_{\mu}$ on
top of the background value $\bar{A}_{\mu}$ satisfying $|\delta
A_{\mu}|\ll|\bar{A}_{\mu}|$ are given by Kase  _et al._ (2018b, c); Kase and
Tsujikawa (2022)

$\delta{A}_{t}=\delta{A}_{r}=0\,,\qquad\delta{A}_{a}=\sum_{l,m}\delta{A}_{lm}(t,r)E_{ab}\,\partial^{b}Y_{lm}(\theta,\varphi)\,,$
(20)

where $\delta{A}_{lm}$ depends on $t$ and $r$.

Under a infinitesimal gauge transformation $x_{\mu}\to x_{\mu}+\xi_{\mu}$,
where $\xi_{t}=0=\xi_{r}$, and
$\xi_{a}=\sum_{l,m}\Lambda_{lm}(t,r)E_{ab}\nabla^{b}Y_{lm}(\theta,\varphi)$,
the metric perturbations transform as $Q_{lm}\to Q_{lm}+\dot{\Lambda}_{lm}$,
$W_{lm}\to W_{lm}+\Lambda_{lm}^{\prime}-2\Lambda_{lm}/r$, and $U_{lm}\to
U_{lm}+2\Lambda_{lm}$ where a dot on top represents the derivative with
respect to $t$. We choose the Regge-Wheeler gauge Regge and Wheeler (1957)
satisfying $\Lambda_{lm}=-U_{lm}/2$ in which $h_{ab}$ in the odd-parity sector
merely vanishes. In the following, we omit the labels $l$ and $m$ of the
quantities $Q_{lm}$, $W_{lm}$, and $\delta{A}_{lm}$ for simplicity.

We expand the action (1) up to quadratic order in odd-parity perturbations and
perform the integration with respect to $\theta$ and $\varphi$. In doing so,
we can focus on the perturbations of $m=0$ mode without loss of generality by
virtue of the underlying background symmetry de Rham  _et al._ (2020). After
repetitively using the integration by parts and the integral formulas for
spherical harmonics (see, e.g., Ref. Zhang  _et al._ (2023a) and references
therein), the second-order action reduces to

$\displaystyle{\cal S}_{\rm odd}^{(2)}$ $\displaystyle=$
$\displaystyle\sum_{l,m}L\int{\rm d}t{\rm d}r\,{\cal L}_{\rm odd}^{(2)}\,,$
(21)

where

$L\equiv l(l+1)\,,$ (22)

and

$\displaystyle{\cal L}_{\rm odd}^{(2)}$ $\displaystyle=$
$\displaystyle\alpha_{1}\left(\dot{W}-Q^{\prime}+\frac{2}{r}Q\right)^{2}+2\left(\alpha_{2}\delta{A}^{\prime}+\alpha_{3}\delta{A}\right)\left(\dot{W}-Q^{\prime}+\frac{2}{r}Q\right)$
(23)
$\displaystyle+\alpha_{4}\dot{\delta{A}}^{2}+\alpha_{5}{\delta{A}}^{\prime
2}+(L-2)\left(\alpha_{6}W^{2}+\alpha_{7}Q^{2}+\alpha_{8}Q\delta{A}\right)+L\alpha_{9}\delta{A}^{2}\,.$

The coefficients $\alpha_{i}$ ($i=1,...,9$) are given by

$\displaystyle\alpha_{1}\equiv\frac{M_{\rm
pl}^{2}\sqrt{h}}{4\sqrt{f}}\,,\qquad\alpha_{2}\equiv\frac{h^{3/2}A_{0}^{\prime}}{2r\sqrt{f}}\left[r\phi^{\prime}f_{3}-4f_{4}+h\phi^{\prime
2}(f_{4,X}+2\tilde{f}_{4})\right]\,,$
$\displaystyle\alpha_{3}\equiv-\frac{\sqrt{h}A_{0}^{\prime}}{2r^{2}\sqrt{f}}\left[r^{2}(f_{2,F}-2h\phi^{\prime
2}f_{2,Y})+4h(r\phi^{\prime}f_{3}-2f_{4})+2h^{2}\phi^{\prime
2}(f_{4,X}+2\tilde{f}_{4})\right]\,,$
$\displaystyle\alpha_{4}\equiv\frac{1}{2r\sqrt{fh}}\left[rf_{2,F}+2h\phi^{\prime}(f_{3}+h\phi^{\prime
2}\tilde{f}_{3})-4h^{\prime}f_{4}+(2h\phi^{\prime\prime}+\phi^{\prime}h^{\prime})\\{rf_{3}+h\phi^{\prime}(f_{4,X}+2\tilde{f}_{4})\\}\right]\,,$
$\displaystyle\alpha_{5}\equiv-\frac{\sqrt{fh}}{2r}\left[r(f_{2,F}-2h\phi^{\prime
2}f_{2,Y})+2h\phi^{\prime}f_{3}+\frac{hf^{\prime}}{f}\\{r\phi^{\prime}f_{3}-4f_{4}+h\phi^{\prime
2}(f_{4,X}+2\tilde{f}_{4})\\}\right]\,,$
$\displaystyle\alpha_{6}\equiv-\frac{\sqrt{fh}}{4r^{2}}\left[M_{\rm
pl}^{2}-\frac{hA_{0}^{\prime 2}}{f}\\{4f_{4}-h\phi^{\prime
2}(f_{4,X}+2\tilde{f}_{4})\\}\right]\,,$
$\displaystyle\alpha_{7}\equiv\frac{1}{4r^{2}\sqrt{fh}}\left(M_{\rm
pl}^{2}-\frac{4hA_{0}^{\prime
2}f_{4}}{f}\right)\,,\qquad\alpha_{8}\equiv\frac{4}{r^{2}}\left(\sqrt{\frac{h}{f}}A_{0}^{\prime}f_{4}\right)^{\prime}\,,$
$\displaystyle\alpha_{9}\equiv-\frac{1}{2r^{2}\sqrt{fh}}\biggl{[}ff_{2,F}-hA_{0}^{\prime
2}f_{2,\tilde{F}\tilde{F}}+2\sqrt{fh}(\sqrt{fh}\phi^{\prime})^{\prime}f_{3}+h^{2}f^{\prime}\phi^{\prime
3}\tilde{f}_{3}-4\sqrt{fh}\left(\sqrt{\frac{h}{f}}f^{\prime}\right)^{\prime}f_{4}$
$\displaystyle\hskip
71.13188pt+h^{3/2}f^{\prime}\phi^{\prime}(\sqrt{h}\phi^{\prime})^{\prime}(f_{4,X}+2\tilde{f}_{4})\biggr{]}\,.$
(24)

Notice that, in the absence of the dependence of $\tilde{F}$ and $Y$ in
$f_{2}$, the above coefficients $\alpha_{i}$ ($i=1,...,9$) reduce to those
defined in Ref. Kase  _et al._ (2018b) up to an overall factor, which are
denoted as $\alpha^{\rm(pre)}_{i}$ ($i=1,...,9$), such that

$\alpha_{i}=r^{2}\sqrt{\frac{f}{h}}\,\alpha^{\rm(pre)}_{i}\,.$ (25)

The second-order action (21) identically vanishes for the monopole mode $l=0$,
i.e., $L=0$. For the dipole mode characterized by $l=1$ ($L=2$), the vector
field perturbation $\delta{A}$ is the only propagating DOF, and it is shown
that the stability conditions for $\delta{A}$ are the same as those for $l\geq
2$ in Ref. Heisenberg  _et al._ (2018). Thus, we focus only on the modes
characterized by $l\geq 2$ and revisit their stability conditions in this
paper.

The dynamical variables in the Lagrangian (23) are $W$ and $\delta{A}$. In
order to integrate out the remaining non-dynamical variable $Q$, we introduce
an auxiliary field $\chi$ to the Lagrangian as

$\displaystyle{\cal L}_{\rm odd}^{(2)}$ $\displaystyle=$
$\displaystyle\alpha_{1}\left[2\chi\left(\dot{W}-Q^{\prime}+\frac{2}{r}Q+\frac{\alpha_{2}\delta{A}^{\prime}+\alpha_{3}\delta{A}}{\alpha_{1}}\right)-\chi^{2}\right]-\frac{(\alpha_{2}\delta{A}^{\prime}+\alpha_{3}\delta{A})^{2}}{\alpha_{1}}$
(26)
$\displaystyle+\alpha_{4}\dot{\delta{A}}^{2}+\alpha_{5}{\delta{A}}^{\prime
2}+(L-2)\left(\alpha_{6}W^{2}+\alpha_{7}Q^{2}+\alpha_{8}Q\delta{A}\right)+L\alpha_{9}\delta{A}^{2}\,.$

Varying the above Lagrangian with respect to $\chi$, we find that the
auxiliary field $\chi$ satisfies

$\chi=\dot{W}-Q^{\prime}+\frac{2}{r}Q+\frac{\alpha_{2}\delta{A}^{\prime}+\alpha_{3}\delta{A}}{\alpha_{1}}\,,$
(27)

which guarantees the equivalence between Eqs. (23) and (26). Meanwhile, the
variation of the Lagrangian (26) with respect to $W$ and $Q$ leads

$\displaystyle\alpha_{1}\dot{\chi}-(L-2)\alpha_{6}W=0\,,$ (28)
$\displaystyle\alpha_{1}\chi^{\prime}+\frac{1}{r^{2}}(r^{2}\alpha_{1})^{\prime}\chi+(L-2)\left(\alpha_{7}Q+\frac{\alpha_{8}}{2}\delta{A}\right)=0\,,$
(29)

respectively. The above algebraic equations can be solved for $W$ and $Q$,
respectively. We substitute these solutions into the Lagrangian (26) and
eliminate the perturbations $Q$ and $W$ from the second-order action. In doing
so, the role of gravitational DOF carried by $Q$ and $W$ in the original
Lagrangian (23) is transferred to the auxiliary field $\chi$ in the resultant
Lagrangian of the form

$(L-2){\cal L}_{\rm
odd}^{(2)}=\dot{\vec{\mathcal{X}}}^{t}{\bm{K}}\dot{\vec{\mathcal{X}}}+\vec{\mathcal{X}}^{\prime
t}{\bm{G}}\vec{\mathcal{X}}^{\prime}+\vec{\mathcal{X}}^{\prime
t}{\bm{S}}\vec{\mathcal{X}}+\vec{\mathcal{X}}^{t}{\bm{M}}\vec{\mathcal{X}}\,,$
(30)

where ${\bm{K,G,S,M}}$ are $2\times 2$ matrices with the non-vanishing
components given by

$\displaystyle K_{11}=-\frac{\alpha_{1}^{2}}{\alpha_{6}}\,,\qquad
K_{22}=(L-2)\alpha_{4}\,,$ $\displaystyle
G_{11}=-\frac{\alpha_{1}^{2}}{\alpha_{7}}\,,\qquad
G_{22}=\frac{(L-2)(\alpha_{1}\alpha_{5}-\alpha_{2}^{2})}{\alpha_{1}}\,,$
$\displaystyle
S_{12}=-S_{21}=-(L-2)\left(\alpha_{2}+\frac{\alpha_{1}\alpha_{8}}{2\alpha_{7}}\right)\,,$
$\displaystyle
M_{11}=-(L-2)\alpha_{1}-\frac{[(r^{2}\alpha_{1})^{\prime}]^{2}}{r^{4}\alpha_{7}}+\left[\frac{\alpha_{1}(r^{2}\alpha_{1})^{\prime}}{r^{2}\alpha_{7}}\right]^{\prime}\,,$
$\displaystyle
M_{22}=-(L-2)\left[\frac{(L-2)\alpha_{8}^{2}}{4\alpha_{7}}-L\alpha_{9}+\frac{\alpha_{3}^{2}}{\alpha_{1}}-\left(\frac{\alpha_{2}\alpha_{3}}{\alpha_{1}}\right)^{\prime}\,\right]\,,$
$\displaystyle
M_{12}=M_{21}=(L-2)\left[\alpha_{3}-\frac{\alpha_{8}(r^{2}\alpha_{1})^{\prime}}{2r^{2}\alpha_{7}}-\frac{1}{2}\left(\alpha_{2}-\frac{\alpha_{1}\alpha_{8}}{2\alpha_{7}}\right)^{\prime}\,\right]\,,$
(31)

and $\vec{\mathcal{X}}$ is the vector defined by

$\vec{\mathcal{X}}^{t}\equiv\left(\chi,\delta{A}\right)\,.$ (32)

The no-ghost conditions, $K_{11}>0$ and $K_{22}>0$, are guaranteed for

$\alpha_{6}<0\,,\qquad\alpha_{4}>0\,.$ (33)

By assuming the solution of the form $\vec{\mathcal{X}}^{t}\propto e^{i(\omega
t-kr)}$ in the small-scale limit characterized by $k\to 0$, the propagation
speed of perturbations along the radial direction is expressed as
$c_{r}=\omega/(\sqrt{fh}\,k)$ in proper time. Substitution of this expression
into the dispersion relation, ${\rm
det}\left(\omega^{2}{\bm{K}}+k^{2}{\bm{G}}\right)=0$, leads

$\displaystyle c_{r1{\rm,odd}}^{2}$ $\displaystyle=$
$\displaystyle-\frac{G_{11}}{fhK_{11}}=-\frac{\alpha_{6}}{fh\,\alpha_{7}}\,,$
(34) $\displaystyle c_{r2{\rm,odd}}^{2}$ $\displaystyle=$
$\displaystyle-\frac{G_{22}}{fhK_{22}}=\frac{\alpha_{2}^{2}-\alpha_{1}\alpha_{5}}{fh\,\alpha_{1}\alpha_{4}}\,,$
(35)

where $c_{r1{\rm,odd}}$ and $c_{r2{\rm,odd}}$ are the propagation speeds of
the odd-parity perturbations arising from the gravity sector and the vector
field perturbation, respectively. On the other hand, by assuming the solution
of the form $\vec{\mathcal{X}}^{t}\propto e^{i(\omega t-l\theta)}$ in the
limit that $L=l(l+1)\gg 1$, the squared propagation speed along the angular
direction is expressed as $c_{\Omega}^{2}=r^{2}\omega^{2}/(l^{2}f)$ in proper
time. We substitute this expression into the dispersion relation of the form
${\rm det}(\omega^{2}{\bm{K}}+{\bm{M}})=0$. Picking up the dominant
contributions for $l\gg 1$, we obtain

$\displaystyle c_{\Omega 1{\rm,odd}}^{2}$ $\displaystyle=$
$\displaystyle-\frac{r^{2}M_{11}}{l^{2}fK_{11}}=-\frac{r^{2}\alpha_{6}}{f\alpha_{1}}\,,$
(36) $\displaystyle c_{\Omega 2{\rm,odd}}^{2}$ $\displaystyle=$
$\displaystyle-\frac{r^{2}M_{22}}{l^{2}fK_{22}}=\frac{r^{2}(\alpha_{8}^{2}-4\alpha_{7}\alpha_{9})}{4f\alpha_{4}\alpha_{7}}\,.$
(37)

Here, $c_{\Omega 1{\rm,odd}}$ is the propagation speed along the angular
direction for the perturbation arising from gravity sector, while $c_{\Omega
2{\rm,odd}}$ is that arising from the vector perturbation. We note that the
expressions of Eqs. (33)-(37) are all the same as those derived in Ref.
Heisenberg  _et al._ (2018) while the dependence of $\tilde{F}$ and $Y$ in
$f_{2}$ alters the coefficients $\alpha_{3}$, $\alpha_{4}$, $\alpha_{5}$, and
$\alpha_{9}$.

For the absence of Laplacian instabilities, we require the conditions
$c_{r1{\rm,odd}}^{2}\geq 0$, $c_{r2{\rm,odd}}^{2}\geq 0$, $c_{\Omega
1{\rm,odd}}^{2}\geq 0$, and $c_{\Omega 2{\rm,odd}}^{2}\geq 0$. Among these
conditions, Eq. (36) shows that the condition $c_{\Omega 1{\rm,odd}}^{2}\geq
0$ is automatically satisfied under the no-ghost condition (33) since
$\alpha_{1}$ is non-negative by definition [cf., (24)]. From Eqs. (34), (35),
(37), the other three conditions for the absence of Laplacian instabilities
are translated to give

$\displaystyle\alpha_{7}\geq 0\,,\qquad\alpha_{2}^{2}-\alpha_{1}\alpha_{5}\geq
0\,,\qquad\alpha_{8}^{2}-4\alpha_{7}\alpha_{9}\geq 0\,.$ (38)

The above stability conditions will assist us below in understanding those for
the even-parity sector.

##  IV stability conditions against even-parity perturbations

We proceed to derive the stability conditions against even-parity
perturbations. On top of the background spacetime characterized by Eq. (8),
the components of metric perturbation $h_{\mu\nu}$ in the even-parity sector
are given by

$\displaystyle
h_{tt}=f(r)\sum_{l,m}H_{0,lm}(t,r)Y_{lm}(\theta,\varphi)\,,\qquad
h_{tr}=h_{rt}=\sum_{l,m}H_{1,lm}(t,r)Y_{lm}(\theta,\varphi)\,,$ $\displaystyle
h_{rr}=h(r)^{-1}\,\sum_{l,m}H_{2,lm}(t,r)Y_{lm}(\theta,\varphi)\,,$
$\displaystyle
h_{ta}=h_{at}=\sum_{l,m}h_{0,lm}(t,r)\nabla_{a}Y_{lm}(\theta,\varphi)\,,\qquad
h_{ra}=h_{ar}=\sum_{l,m}h_{1,lm}(t,r)\nabla_{a}Y_{lm}(\theta,\varphi)\,,$
$\displaystyle
h_{ab}=\sum_{l,m}\left[K_{lm}(t,r)g_{ab}Y_{lm}(\theta,\varphi)+G_{lm}(t,r)\nabla_{a}\nabla_{b}Y_{lm}(\theta,\varphi)\right]\,,$
(39)

where $H_{0,lm}$, $H_{1,lm}$, $H_{2,lm}$, $h_{0,lm}$, $h_{1,lm}$, $K_{lm}$,
and $G_{lm}$ are scalar quantities depending on $t$ and $r$ [One should avoid
confusing the $K$ and $G$ in here with those appearing in, e.g., (31)].
Similarly, we consider the perturbations of scalar and vector fields on top of
their background values as444 We note that we shall in general omit the
overbar in writing background quantities like $\bar{\phi}$ and
${\bar{A}}_{\mu}$ in the following, since the order of a quantity can be
easily identified by the context and counting the perturbation terms, just
like what we have done in, e.g., Eq.(12).

$\displaystyle\phi$ $\displaystyle=$
$\displaystyle\bar{\phi}+\sum_{l,m}\delta\phi_{lm}(t,r)Y_{lm}(\theta,\varphi)\,,$
(40) $\displaystyle A_{\mu}$ $\displaystyle=$
$\displaystyle\bar{A}_{\mu}+\delta A_{\mu}\,,$ (41)

with

$\delta A_{t}=\sum_{l,m}\delta
A_{0,lm}(t,r)Y_{lm}(\theta,\varphi)\,,\qquad\delta A_{r}=\sum_{l,m}\delta
A_{1,lm}(t,r)Y_{lm}(\theta,\varphi)\,,\qquad\delta A_{a}=\sum_{l,m}\delta
A_{2,lm}(t,r)\nabla_{a}Y_{lm}(\theta,\varphi)\,,$ (42)

where $\bar{\phi}$ and $\bar{A}_{\mu}$ are background values given by Eqs. (9)
and (10), respectively. The scalar quantities $\delta\phi_{lm}$, $\delta
A_{0,lm}$, $\delta A_{1,lm}$, and $\delta A_{2,lm}$ are functions of $t$ and
$r$, which are, as usual, assumed to be much smaller than the background
quantities.

Let us consider the infinitesimal gauge transformation $x_{\mu}\to
x_{\mu}+\xi_{\mu}$ with

$\xi_{t}=\sum_{l,m}{\cal
T}_{lm}(t,r)Y_{lm}(\theta,\varphi)\,,\qquad\xi_{r}=\sum_{l,m}{\cal
R}_{lm}(t,r)Y_{lm}(\theta,\varphi)\,,\qquad\xi_{a}=\sum_{l,m}\Theta_{lm}(t,r)\nabla_{a}Y_{lm}(\theta,\varphi)\,,$
(43)

where ${\cal T}_{lm}$, ${\cal R}_{lm}$, and $\Theta_{lm}$ are the scalar
quantities depending on $t$ and $r$. In the following we omit the subscripts
$l$ and $m$ of the scalar quantities in Eqs. (39) and (43) for simplicity.
Under the above transformation, the scalar quantities in Eq. (39) transform as

$\displaystyle H_{0}\to H_{0}+\frac{2}{f}\dot{{\cal
T}}-\frac{f^{\prime}h}{f}{\cal R}\,,\qquad H_{1}\to H_{1}+\dot{{\cal R}}+{\cal
T}^{\prime}-\frac{f^{\prime}}{f}{\cal T}\,,\qquad H_{2}\to H_{2}+2h{\cal
R}^{\prime}+h^{\prime}{\cal R}\,,$ $\displaystyle h_{0}\to h_{0}+{\cal
T}+\dot{\Theta}\,,\qquad h_{1}\to h_{1}+{\cal
R}+\Theta^{\prime}-\frac{2}{r}\Theta\,,\qquad K\to K+\frac{2}{r}h{\cal
R}\,,\qquad G\to G+\frac{2}{r^{2}}\Theta\,,$
$\displaystyle\delta\phi\to\delta\phi-\phi^{\prime}h{\cal R}\,,$ (44)

where we have dropped the $lm$ in subscripts for the corresponding quantities
for simplicity. For $l\geq 2$, we choose the gauge given by ${\cal
T}=-h_{0}+{r^{2}}\dot{G}/2$, ${\cal R}=-rK/({2h})$, and $\Theta=-{r^{2}}G/2$
so that the perturbations $h_{0}$, $K$, $G$ 555One should avoid confusing the
$G$ in here with the gravitational constant mentioned previously.  identically
vanish. This gauge fixing is equivalent to setting

$h_{0}=0\,,\qquad K=0\,,\qquad G=0\,,$ (45)

in Eq. (39) from the beginning (This is sometimes referred as the EZ gauge
Thompson  _et al._ (2017)). We also consider the $U(1)$ gauge transformation,

$\delta A_{\mu}\to\delta A_{\mu}+\partial_{\mu}\delta\chi\qquad{\rm
with}\qquad\delta\chi=\sum_{l,m}\tilde{\chi}(t,r)Y_{lm}(\theta,\varphi)\,,$
(46)

under which the scalar quantities in Eq. (42) transform as

$\delta A_{0}\to\delta A_{0}+\dot{\tilde{\chi}}\,,\qquad\delta A_{1}\to\delta
A_{1}+\tilde{\chi}^{\prime}\,,\qquad\delta A_{2}\to\delta
A_{2}+\tilde{\chi}\,.$ (47)

For the gauge choice $\tilde{\chi}=-\delta A_{2}$, the quantity $\delta A_{2}$
identically vanishes. This corresponds to setting

$\delta A_{2}=0\,,$ (48)

in Eq. (42) from the beginning. Thus, with all the gauge choices mentioned
above, we are now left with 7 variables, i.e., $\\{\delta\phi,\delta
A_{0},\delta A_{1},h_{1},H_{0},H_{1},H_{2}\\}$ (and at this point they
themselves, as well as their combinations, represent gauge invariants).

###  IV.1 Second-order action and perturbation equations of motion

We expand the action (1) up to second order in terms of the even-parity
perturbations given in Eqs. (39), (40), and (41) under the gauge choices (45)
and (48). In doing so, we can focus on the $m=0$ mode without loss of
generality for the same reason as in the case of odd-parity perturbations.
After lengthy but straightforward calculation, the second-order action in the
even-parity sector reduces to

${\cal S}_{\rm even}=\sum_{l}\int{\rm d}t\,{\rm d}r\left({\cal L}_{u}+{\cal
L}_{A}\right)\,,$ (49)

where

$\displaystyle{\cal L}_{u}$ $\displaystyle\equiv$ $\displaystyle
H_{0}\left[a_{1}\delta\phi^{\prime\prime}+a_{2}\delta\phi^{\prime}+a_{3}H_{2}^{\prime}+La_{4}h_{1}^{\prime}+\left(a_{5}+La_{6}\right)\delta\phi+\left(a_{7}+La_{8}\right)H_{2}+La_{9}h_{1}\right]$
(50)
$\displaystyle+Lb_{1}H_{1}^{2}+H_{1}\left(b_{2}\dot{\delta\phi}^{\prime}+b_{3}\dot{\delta\phi}+b_{4}\dot{H}_{2}+Lb_{5}\dot{h}_{1}\right)+c_{1}\dot{\delta\phi}\dot{H}_{2}+H_{2}\left[c_{2}\delta\phi^{\prime}+(c_{3}+Lc_{4})\delta\phi+Lc_{5}h_{1}\right]$
$\displaystyle+c_{6}H_{2}^{2}+Ld_{1}\dot{h}_{1}^{2}+Lh_{1}\left(d_{2}\delta\phi^{\prime}+d_{3}\delta\phi\right)+Ld_{4}h_{1}^{2}+e_{1}\dot{\delta\phi}^{2}+e_{2}\delta\phi^{\prime
2}+\left(e_{3}+Le_{4}\right)\delta\phi^{2}\,,$ $\displaystyle{\cal L}_{A}$
$\displaystyle\equiv$ $\displaystyle v_{1}(\delta A_{0}^{\prime}-\dot{\delta
A}_{1})^{2}+(\delta A_{0}^{\prime}-\dot{\delta
A}_{1})\left(v_{2}H_{0}+v_{3}H_{2}+v_{4}\delta\phi^{\prime}+v_{5}\delta\phi+Lv_{6}h_{1}\right)+\frac{1}{2}Lv_{6}h_{1}\dot{\delta
A}_{1}+v_{7}H_{0}^{2}$ (51) $\displaystyle+L(v_{8}h_{1}\delta
A_{0}+v_{9}\delta A_{0}^{2}+v_{10}\delta A_{1}^{2}+v_{11}H_{1}\delta
A_{1}+v_{12}H_{2}\delta A_{0}+v_{13}\delta\phi\delta A_{0})\,.$

The coefficients $a_{1}$, $a_{2}$, …, $v_{9}$ are given in Appendix A. In Ref.
Kase and Tsujikawa (2023), the even-parity perturbations in Horndeski theories
with the interaction between scalar and Maxwell fields characterized by the
Lagrangian $G_{2}(\phi,X,F)$ are studied. Compared to that, the presence of
the $U(1)$ GI SVT interactions characterized by $f_{3}$, $\tilde{f}_{3}$,
$f_{4}$, and $\tilde{f}_{4}$ in Eq. (1) gives rise to the new terms with the
coefficients $v_{6}$, $v_{11}$, $v_{12}$, and $v_{13}$, in Eq. (51). For
$f_{2}=G_{2}(\phi,X,F)$, $f_{3}=0$, $\tilde{f}_{3}=0$, $f_{4}=0$, and
$\tilde{f}_{4}=0$, these coefficients identically vanish and the Lagrangian
${\cal L}_{A}$ coincides with that derived in Ref. Gannouji and Baez (2022);
Kase and Tsujikawa (2023). Since the Lagrangian ${\cal L}_{A}$ possesses the
similar structure to that in Ref. Kase and Tsujikawa (2023), we can resort to
the analogous method in order to identify the dynamical vector DOF in the
even-parity sector. Let us introduce an auxiliary field $V(t,r)$ and rewrite
Eq. (51) as follows,

$\displaystyle{\cal L}_{A}$ $\displaystyle=$ $\displaystyle
v_{1}\left\\{2V\left[\delta A_{0}^{\prime}-\dot{\delta
A}_{1}+\frac{1}{2v_{1}}\left(v_{2}H_{0}+v_{3}H_{2}+v_{4}\delta\phi^{\prime}+v_{5}\delta\phi+Lv_{6}h_{1}\right)\right]-V^{2}\right\\}$
(52)
$\displaystyle-\frac{1}{4v_{1}}\left(v_{2}H_{0}+v_{3}H_{2}+v_{4}\delta\phi^{\prime}+v_{5}\delta\phi+Lv_{6}h_{1}\right)^{2}+\frac{1}{2}Lv_{6}h_{1}\dot{\delta
A}_{1}+v_{7}H_{0}^{2}$ $\displaystyle+L(v_{8}h_{1}\delta A_{0}+v_{9}\delta
A_{0}^{2}+v_{10}\delta A_{1}^{2}+v_{11}H_{1}\delta A_{1}+v_{12}H_{2}\delta
A_{0}+v_{13}\delta\phi\delta A_{0})\,.$

By varying this action with respect to $V$, we find that

$\displaystyle V=\delta A_{0}^{\prime}-\dot{\delta
A}_{1}+\frac{1}{2v_{1}}\left(v_{2}H_{0}+v_{3}H_{2}+v_{4}\delta\phi^{\prime}+v_{5}\delta\phi+Lv_{6}h_{1}\right)\,,$
(53)

which shows the equivalence between Eqs. (51) and (52). Nevertheless, the
auxiliary field $V$ plays a key role to find out the dynamical vector DOF on
the point that the quadratic derivative terms $\delta A_{0}^{\prime 2}$ and
$\dot{\delta A}_{1}^{2}$ in Eq. (51) do not appear in Eq. (52). The
disappearance of these derivative terms in Eq. (52) allows us to solve the
perturbation equations of $\delta A_{0}$ and $\delta A_{1}$ for themselves
explicitly in analogy with what we did in Eqs. (28)-(29). On using these
solutions, the dynamical property of vector field perturbations can be
aggregated into the auxiliary field $V$ as we will see later. Moreover, the
quadratic term $H_{0}^{2}$ in Eq. (52) identically vanishes by virtue of the
relation among the coefficients $v_{1}$, $v_{2}$, and $v_{7}$, of the form,

$v_{7}=\frac{v_{2}^{2}}{4v_{1}}\,.$ (54)

This shows that the perturbation $H_{0}$ appears only linearly in the total
action (49) with Eqs. (50) and (52). In other words, the perturbation $H_{0}$
corresponds to a Lagrange multiplier. The variation of the action with respect
to $H_{0}$ gives constraint on other perturbation variables, and $H_{0}$
simply disappears once the constraint is applied to the action.

We vary the total action (49) represented by Eqs. (50) and (52) with respect
to $H_{0}$, $H_{1}$, $H_{2}$, $h_{1}$, $\delta A_{0}$, $\delta A_{1}$, and
$\delta\phi$, so as to obtain the following linear perturbation equations,

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
a_{3}H_{2}^{\prime}+La_{4}h_{1}^{\prime}+\left(a_{2}-\frac{v_{2}v_{4}}{2v_{1}}\right)\delta\phi^{\prime}+\left(La_{6}+a_{5}-\frac{v_{2}v_{5}}{2v_{1}}\right)\delta\phi+\left(La_{8}+a_{7}-\frac{v_{2}v_{3}}{2v_{1}}\right)H_{2}$
(55)
$\displaystyle+L\left(a_{9}-\frac{v_{2}v_{6}}{2v_{1}}\right)h_{1}+v_{2}V\,,$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
2Lb_{1}H_{1}+b_{3}\dot{\delta\phi}+b_{4}\dot{H_{2}}+Lb_{5}\dot{h_{1}}+Lv_{11}\delta
A_{1}\,,$ (56) $\displaystyle 0$ $\displaystyle=$ $\displaystyle-
b_{4}\dot{H}_{1}+\left(c_{2}-\frac{v_{3}v_{4}}{2v_{1}}\right)\delta\phi^{\prime}+\left(c_{3}+Lc_{4}-\frac{v_{3}v_{5}}{2v_{1}}\right)\delta\phi+L\left(c_{5}-\frac{v_{3}v_{6}}{2v_{1}}\right)h_{1}+\left(2c_{6}-\frac{v_{3}^{2}}{2v_{1}}\right)H_{2}$
(57) $\displaystyle-
a_{3}H_{0}^{\prime}+\left(a_{7}-a_{3}^{\prime}+La_{8}-\frac{v_{2}v_{3}}{2v_{1}}\right)H_{0}+v_{3}V+Lv_{12}\delta
A_{0}\,,$ $\displaystyle 0$ $\displaystyle=$
$\displaystyle-2d_{1}\ddot{h_{1}}+\left(d_{2}-\frac{v_{4}v_{6}}{2v_{1}}\right)\delta\phi^{\prime}+\left(d_{3}-\frac{v_{5}v_{6}}{2v_{1}}\right)\delta\phi+\left(2d_{4}-\frac{Lv_{6}^{2}}{2v_{1}}\right)h_{1}-a_{4}H_{0}^{\prime}+\left(a_{9}-a_{4}^{\prime}-\frac{v_{2}v_{6}}{2v_{1}}\right)H_{0}$
(58) $\displaystyle-
b_{5}\dot{H_{1}}+\left(c_{5}-\frac{v_{3}v_{6}}{2v_{1}}\right)H_{2}+v_{6}V+\frac{1}{2}v_{6}\dot{\delta
A}_{1}+v_{8}\delta A_{0}\,,$ $\displaystyle 0$ $\displaystyle=$
$\displaystyle-2(v_{1}V)^{\prime}+L\left(v_{8}h_{1}+2v_{9}\delta
A_{0}+v_{12}H_{2}+v_{13}\delta\phi\right)\,,$ (59) $\displaystyle 0$
$\displaystyle=$ $\displaystyle
2v_{1}\dot{V}+L\left(-\frac{1}{2}v_{6}\dot{h_{1}}+2v_{10}\delta
A_{1}+v_{11}H_{1}\right)\,,$ (60) $\displaystyle 0$ $\displaystyle=$
$\displaystyle-2e_{1}\ddot{\delta\phi}-\left(2e_{2}-\frac{v_{4}^{2}}{2v_{1}}\right)\delta\phi^{\prime\prime}+\left[2e_{3}+2Le_{4}+\left(\frac{v_{4}v_{5}}{2v_{1}}\right)^{\prime}-\frac{v_{5}^{2}}{2v_{1}}\right]\delta\phi-\left(a_{2}-\frac{v_{2}v_{4}}{2v_{1}}\right)H_{0}^{\prime}$
(61)
$\displaystyle+\left[a_{5}+La_{6}-a_{2}^{\prime}+\left(\frac{v_{2}v_{4}}{2v_{1}}\right)^{\prime}-\frac{v_{2}v_{5}}{2v_{1}}\right]H_{0}-b_{3}\dot{H}_{1}-\left(c_{2}-\frac{v_{3}v_{4}}{2v_{1}}\right)H_{2}^{\prime}$
$\displaystyle+\left[c_{3}+Lc_{4}-c_{2}^{\prime}+\left(\frac{v_{3}v_{4}}{2v_{1}}\right)^{\prime}-\frac{v_{3}v_{5}}{2v_{1}}\right]H_{2}-L\left(d_{2}-\frac{v_{4}v_{6}}{2v_{1}}\right)h_{1}^{\prime}+L\left[d_{3}-d_{2}^{\prime}+\left(\frac{v_{4}v_{6}}{2v_{1}}\right)^{\prime}-\frac{v_{5}v_{6}}{2v_{1}}\right]h_{1}$
$\displaystyle-\left[2e_{2}^{\prime}-\left(\frac{v_{4}^{2}}{2v_{1}}\right)^{\prime}\right]\delta\phi^{\prime}-v_{4}V^{\prime}+(v_{5}-v_{4}^{\prime})V+Lv_{13}\delta
A_{0}\,,$

where we used the relation (54) and substituted in the coefficients being 0 in
Appendix A. In the following, we study the linear stability conditions for the
three cases (1) $l\geq 2$, (2) $l=0$, and (3) $l=1$, in turn.

###  IV.2 Linear stability conditions for $l\geq 2$

We introduce the following quantity

$\displaystyle\psi$ $\displaystyle\equiv$ $\displaystyle
H_{2}-\frac{L}{r}h_{1}\,,$ (62)

which corresponds to the propagating DOF of gravitational sector as in the
case of GR. Now, we have all the representations for the 3 DOFs, i.e.,
gravitational ($\psi$), vector ($V$) and scalar ($\delta\phi$), in hand. We
replace the $H_{2}$, $\dot{H}_{2}$, and $H_{2}^{\prime}$ in (55)-(60) with
$\psi$ and its derivatives. In doing so, the quantity $h_{1}^{\prime}$
disappears from Eq. (55) on account of the relation $a_{3}=-ra_{4}$, and the
resultant equation can be explicitly solved for $h_{1}$ as a combination of
$\psi$, $\delta\phi^{\prime}$, $\delta\phi$, and $V$. Substituting this
solution into Eqs. (56), (59), (60) and combining them, we can also express
$H_{1}$, $\delta A_{0}$, and $\delta A_{1}$ in terms of $\psi$, $\delta\phi$,
$V$, and their derivatives. Finally, by substituting these solutions back into
the total action (49) with Eqs. (50) and (52), the perturbations $H_{1}$,
$H_{2}$, $h_{1}$, $\delta A_{0}$, $\delta A_{1}$ are removed. The quantity
$H_{0}$ is also removed from the action as we discussed below Eq. (54). As a
consequence, the resultant action is composed of the only 3 dynamical
perturbations $\psi$, $\delta\phi$, and $V$. Denoting the reduced Lagrangian
as ${\cal L}_{\rm even}^{(2)}$ by mimicking (30), it can be written as

$\displaystyle{\cal L}_{\rm
even}^{(2)}=\dot{\vec{\mathcal{Y}}}^{t}{\tilde{\bm{K}}}\dot{\vec{\mathcal{Y}}}+\vec{\mathcal{Y}}^{\prime
t}{\tilde{\bm{G}}}\vec{\mathcal{Y}}^{\prime}+\vec{\mathcal{Y}}^{\prime
t}{\tilde{\bm{S}}}\vec{\mathcal{Y}}+\vec{\mathcal{Y}}^{t}{\tilde{\bm{M}}}\vec{\mathcal{Y}}\,,$
(63)

where $\vec{\mathcal{Y}}$ is the vector defined by

$\vec{\mathcal{Y}}^{t}=\left(\psi,\delta\phi,V\right)\,.$ (64)

Here, the matrices ${\tilde{\bm{K}}}$, ${\tilde{\bm{G}}}$ and
${\tilde{\bm{M}}}$ are symmetric while ${\tilde{\bm{S}}}$ is anti-symmetric.
Putting together the 2 DOFs from the odd-parity sector [cf., (32)] and the 3
DOFs from the even-parity sector [cf., (64)], we obtain 5 DOFs in total, as
once promised.

####  IV.2.1 No-ghost conditions

The components of kinetic matrix ${\tilde{\bm{K}}}$ are given as

$\displaystyle\tilde{K}_{11}=\frac{1}{L}\left[{\cal
K}_{1}+\left(\frac{a_{4}^{\prime}}{a_{4}}+\frac{1}{r}+\frac{L}{2rh}\right){\cal
K}_{2}-\frac{1}{2}{\cal K}_{2}^{\prime}\right]\,,$
$\displaystyle\tilde{K}_{12}=\frac{1}{2L}\left[\frac{fb_{3}}{ra_{4}}\left\\{{\cal
K}_{1}+\frac{1}{2}\left(\frac{a_{4}^{\prime}}{a_{4}}+\frac{1}{r}+\frac{L}{2rh}\right){\cal
K}_{2}\right\\}-\frac{1}{ra_{4}}\left(a_{5}+La_{6}-\frac{A_{0}^{\prime}}{2}v_{5}\right){\cal
K}_{2}-\frac{1}{2}\left(\frac{fb_{3}}{ra_{4}}{\cal
K}_{2}\right)^{\prime}\right]\,,$
$\displaystyle\tilde{K}_{22}=e_{1}+\frac{1}{4L}\left[\frac{fb_{3}}{r^{2}a_{4}^{2}}\left\\{fb_{3}{\cal
K}_{1}-2\left(a_{5}+La_{6}-\frac{A_{0}^{\prime}}{2}v_{5}\right){\cal
K}_{2}\right\\}-\frac{1}{2}\left\\{\left(\frac{fb_{3}}{ra_{4}}\right)^{2}{\cal
K}_{2}\right\\}^{\prime}\,\right]\,,$
$\displaystyle\tilde{K}_{13}=\frac{v_{1}}{4Lra_{4}}\left(\frac{fv_{6}}{v_{10}}{\cal
K}_{1}-2A_{0}^{\prime}{\cal
K}_{2}\right)\,,\qquad\tilde{K}_{23}=\frac{fb_{3}}{2ra_{4}}\tilde{K}_{13}\,,\qquad\tilde{K}_{33}=\frac{fv_{1}^{2}}{2Lr^{2}a_{4}v_{10}}{\cal
K}_{1}\,,$ (65)

where we used the relations among the coefficients given in Appendix A, and
introduced

${\cal
K}_{1}\equiv-\frac{2r^{2}a_{4}}{f}\left(1-\frac{fv_{6}^{2}}{8a_{4}v_{10}}\right)^{-1}\,,\qquad{\cal
K}_{2}\equiv\frac{4r^{2}a_{4}}{f}\left(\frac{f^{\prime}}{f}+\frac{L}{rh}-\frac{2}{r}+\frac{A_{0}^{\prime}v_{6}}{2a_{4}}\right)^{-1}\,.$
(66)

The ghost instabilities are absent if the kinetic matrix ${\tilde{\bm{K}}}$ is
positive definite, which is translated to requiring the following three no-
ghost conditions,

$\displaystyle{\tilde{K}}_{33}$ $\displaystyle>$ $\displaystyle 0\,,$ (67)
$\displaystyle{\tilde{K}}_{11}{\tilde{K}}_{33}-{\tilde{K}}_{13}{\tilde{K}}_{31}$
$\displaystyle>$ $\displaystyle 0\,,$ (68)
$\displaystyle\det{{\tilde{\bm{K}}}}$ $\displaystyle>$ $\displaystyle 0\,.$
(69)

On using the relations among the coefficients shown in Appendix A, the
expression of the first no-ghost condition (67) reduces to

${\tilde{K}}_{33}=\frac{\alpha_{1}v_{1}^{2}}{L(\alpha_{2}^{2}-\alpha_{1}\alpha_{5})}>0\,.$
(70)

Remembering that $\alpha_{1}$ is positive by definition, this shows that the
first no-ghost condition automatically gets satisfied provided the odd-parity
stability conditions given in Eq. (38). Although the remaining two no-ghost
conditions are slightly complicated, they can be simplified by the use of the
following relation,

${\cal K}_{2}^{\prime}=\frac{(r^{2}fh)^{\prime}}{r^{2}fh}{\cal
K}_{2}-\frac{1}{8rf^{2}a_{4}}\left[L\left(\frac{f^{3}}{r^{2}h}\right)^{\prime}+\frac{2f(rf^{\prime}-2f)^{2}}{r^{3}}-\frac{4f^{2}(2f^{2}\alpha_{7}+A_{0}^{\prime
2}v_{10})}{ra_{4}}\right]{\cal K}_{2}^{2}\,,$ (71)

where we used Eq. (131) in Appendix A for the derivation of this relation.
Eliminating ${\cal K}_{2}^{\prime}$ in Eq. (65) by using the above relation
and substituting the definition of ${\cal K}_{1}$ and ${\cal K}_{2}$, the
second and third no-ghost conditions (68)-(69) reduce to

$\displaystyle{\tilde{K}}_{11}{\tilde{K}}_{33}-{\tilde{K}}_{13}{\tilde{K}}_{31}=\frac{8r^{4}hv_{1}^{2}a_{4}^{3}[L(2\sqrt{fh}{\cal
P}_{1}+rA_{0}^{\prime}v_{6})-8r^{2}fh\alpha_{7}]}{L^{2}f^{2}(4Lra_{4}+2r^{2}hA_{0}^{\prime}v_{6}+\sqrt{fh}{\cal
P}_{2})^{2}(\alpha_{2}^{2}-\alpha_{1}\alpha_{5})}$ $\displaystyle>$
$\displaystyle 0\,,$ (72)
$\displaystyle\det{\tilde{\bm{K}}}=\frac{32r^{6}hv_{1}^{2}a_{4}^{3}\alpha_{7}(L-2)(2\sqrt{fh}{\cal
P}_{1}+rA_{0}^{\prime}v_{6}-4r^{2}fh\alpha_{7})}{L^{2}f^{2}\phi^{\prime
2}(4Lra_{4}+2r^{2}hA_{0}^{\prime}v_{6}+\sqrt{fh}{\cal
P}_{2})^{2}(\alpha_{2}^{2}-\alpha_{1}\alpha_{5})}$ $\displaystyle>$
$\displaystyle 0\,,$ (73)

respectively, where

$\displaystyle{\cal
P}_{1}\equiv\frac{1}{f}\left(\frac{r\sqrt{f}}{\sqrt{h}}\right)^{\prime}a_{4}\,,\qquad{\cal
P}_{2}\equiv\frac{8r^{3}\sqrt{h}}{f}\left(\frac{\sqrt{f}}{r}\right)^{\prime}a_{4}\,.$
(74)

Under the stability conditions in the odd-parity sector given in Eq. (38), the
no-ghost conditions (72) and (73) are satisfied if $L(2\sqrt{fh}{\cal
P}_{1}+rA_{0}^{\prime}v_{6})-8r^{2}fh\alpha_{7}>0$ and $2\sqrt{fh}{\cal
P}_{1}+rA_{0}^{\prime}v_{6}-4r^{2}fh\alpha_{7}>0$ hold, respectively. Since
$L\geq 2$ and $\alpha_{7}\geq 0$, we notice that the former is always ensured
as long as the latter holds. Hence, the second and third no-ghost conditions
converge to

$\mathcal{K}\equiv 2\left({\cal
P}_{1}+\frac{v_{6}A_{0}^{\prime}r}{2\sqrt{fh}}\right)-4\sqrt{fh}\alpha_{7}r^{2}>0\,.$
(75)

In conclusion, the no-ghost conditions in the even-parity sector give rise to
only one additional constraint (75) to the stability conditions in the odd-
parity sector.

####  IV.2.2 Radial and angular Laplacian stability conditions

By assuming the solution of the form $\vec{\mathcal{Y}}^{t}\propto e^{i(\omega
t-kr)}$ in the small-scale limit characterized by $k\to 0$, the propagation
speed of perturbations along the radial direction is expressed as
$c_{r}=\omega/(\sqrt{fh}\,k)$ in proper time. Substitution of this expression
into the dispersion relation reads

${\rm det}\left(fhc_{r}^{2}{\tilde{\bm{K}}}+{\tilde{\bm{G}}}\right)=0\,.$ (76)

The components of matrix ${\tilde{\bm{G}}}$ are given as

$\displaystyle\tilde{G}_{11}=-\frac{f^{2}{\cal
K}_{2}^{2}}{16r^{4}a_{4}^{2}}\left[\frac{1}{v_{1}}\left(A_{0}^{\prime}v_{1}+\frac{\phi^{\prime}}{2}v_{4}\right)^{2}-4(rc_{5}+c_{6})-\frac{4r^{2}}{L}d_{4}+\frac{(2rhv_{8}+Lv_{6})^{2}}{4Lh^{2}v_{9}}\right]\,,$
$\displaystyle\tilde{G}_{12}=\frac{A_{0}^{\prime}v_{4}-2a_{2}}{2ra_{4}}\tilde{G}_{11}-\frac{f{\cal
K}_{2}}{4r^{4}a_{4}}\left[c_{2}+rd_{2}+\frac{v_{4}}{2v_{1}}\left(A_{0}^{\prime}v_{1}+\frac{\phi^{\prime}}{2}v_{4}\right)\right]\,,$
$\displaystyle\tilde{G}_{22}=e_{2}-\frac{v_{4}^{2}}{4v_{1}}+\frac{A_{0}^{\prime}v_{4}-2a_{2}}{ra_{4}}\left(\tilde{G}_{12}-\frac{A_{0}^{\prime}v_{4}-2a_{2}}{4ra_{4}}\tilde{G}_{11}\right)\,,\qquad\frac{\tilde{G}_{13}}{\tilde{K}_{13}}=\frac{\tilde{G}_{23}}{\tilde{K}_{23}}=\frac{\tilde{G}_{33}}{\tilde{K}_{33}}=-\frac{fv_{6}^{2}-8a_{4}v_{10}}{8a_{4}v_{9}}\,.$
(77)

By virtue of the relation between $\tilde{G}_{i3}$ and $\tilde{K}_{i3}$
($i=1,2,3$), the propagation speed of the vector DOF decouples from the other
two in the dispersion relation (76). On using the relations among the
coefficients given in Appendix A, it reduces to

$c_{r3{\rm,even}}^{2}=\frac{\alpha_{2}^{2}-\alpha_{1}\alpha_{5}}{fh\alpha_{1}\alpha_{4}}\,,$
(78)

Although the remaining two propagation speeds look being coupled with each
other in the dispersion relation (76), they decouples on the use of the
relations among coefficients given in Appendix A such that

$\displaystyle c_{r1{\rm,even}}^{2}$ $\displaystyle=$
$\displaystyle-\frac{\alpha_{6}}{fh\alpha_{7}}\,,$ (79) $\displaystyle
c_{r2{\rm,even}}^{2}$ $\displaystyle=$
$\displaystyle\frac{\phi^{\prime}\left[v_{4}\left(2A_{0}^{\prime}v_{1}+\phi^{\prime}v_{4}\right)+4v_{1}(c_{2}+4rd_{2})\right]}{4v_{1}(2\sqrt{fh}{\cal
P}_{1}+rA_{0}^{\prime}v_{6}-4r^{2}fh\alpha_{7})}\,.$ (80)

They correspond to the radial propagation speeds of the tenor and scalar mode,
respectively. For the absence of Laplacian instabilities along the radial
direction, we require

$c_{r1{\rm,even}}^{2}\geq 0\,,\qquad c_{r2{\rm,even}}^{2}\geq 0\,,\qquad
c_{r3{\rm,even}}^{2}\geq 0\,.$ (81)

By observation, we notice that the squared radial propagation speeds of vector
(78) and tensor DOFs (79) are equivalent to those in the odd-parity sector
given by (35) and (34), respectively. Thus, the first and third radial
Laplacian stability conditions are automatically satisfied as long as the odd-
parity sector is stable.

On the other hand, by assuming the solution of the form
$\vec{\mathcal{Y}}^{t}\propto e^{i(\omega t-l\theta)}$ in the limit that
$L=l(l+1)\gg 1$, the squared propagation speed along the angular direction
(which is dimensionless) is expressed as
$c_{\Omega}^{2}=r^{2}\omega^{2}/(l^{2}f)$ in proper time. We substitute this
expression into the dispersion relation of the form

${\rm
det}\left(\frac{fl^{2}c_{\Omega}^{2}}{r^{2}}{\tilde{\bm{K}}}+{\tilde{\bm{M}}}\right)=0\,.$
(82)

We expand this equation for $L\gg 1$ and pick up the dominant contribution so
as to derive the propagation speeds along the angular direction. As we will
see later, the leading order contribution possesses the linear dependence in
$L$, but this contribution identically vanishes after the substitution of
relations among coefficients. This shows that we need to extract the sub-
leading order contribution from Eq. (82). In order to do so, we expand each
component of matrices $\tilde{\bm{K}}$ and $\tilde{\bm{M}}$ for $L\gg 1$ up to
the sub-leading order and find that the components possess the following
$L$-dependence,

$\displaystyle\tilde{K}_{11}=\frac{\tilde{K}^{(0)}_{11}}{L}+\frac{\tilde{K}^{(1)}_{11}}{L^{2}}\,,\qquad\tilde{K}_{12}=\frac{\tilde{K}^{(0)}_{12}}{L}+\frac{\tilde{K}^{(1)}_{12}}{L^{2}}\,,\qquad\tilde{K}_{13}=\frac{\tilde{K}^{(0)}_{13}}{L}+\frac{\tilde{K}^{(1)}_{13}}{L^{2}}\,,$
$\displaystyle\tilde{K}_{22}=\tilde{K}^{(0)}_{22}+\frac{\tilde{K}^{(1)}_{22}}{L}\,,\qquad\tilde{K}_{23}=\frac{\tilde{K}^{(0)}_{23}}{L}+\frac{\tilde{K}^{(1)}_{23}}{L^{2}}\,,\qquad\tilde{K}_{33}=\frac{\tilde{K}^{(0)}_{33}}{L}\,,$
$\displaystyle\tilde{M}_{11}=\tilde{M}^{(0)}_{11}+\frac{\tilde{M}^{(1)}_{11}}{L}\,,\qquad\tilde{M}_{12}=\tilde{M}^{(0)}_{12}+\frac{\tilde{M}^{(1)}_{12}}{L}\,,\qquad\tilde{M}_{13}=\tilde{M}^{(0)}_{13}+\frac{\tilde{M}^{(1)}_{13}}{L}\,,$
$\displaystyle\tilde{M}_{22}=L\tilde{M}^{(0)}_{22}+\tilde{M}^{(1)}_{22}\,,\qquad\tilde{M}_{23}=\tilde{M}^{(0)}_{23}+\frac{\tilde{M}^{(1)}_{23}}{L}\,,\qquad\tilde{M}_{33}=\tilde{M}^{(0)}_{33}+\frac{\tilde{M}^{(1)}_{33}}{L}\,,$
(83)

where the quantities with superscripts $(0)$ and $(1)$ correspond to the
coefficients of leading and sub-leading contributions for $L\gg 1$,
respectively, which do not contain $L$ themselves. We note that
$\tilde{K}_{33}$ cannot be expanded further as seen from Eq. (65).
Substituting these expressions into the dispersion relation (82) and expanding
it, the leading order contribution is in proportion to $L$. As mentioned
above, this contribution identically vanishes since the quantities
$\tilde{M}^{(0)}_{ij}$ $(i,j=1,2,3)$ satisfy

$\frac{\tilde{M}^{(0)}_{11}}{\tilde{K}^{(0)}_{11}}=\frac{\tilde{M}^{(0)}_{13}}{\tilde{K}^{(0)}_{13}}=\frac{\tilde{M}^{(0)}_{33}}{\tilde{K}^{(0)}_{33}}=\frac{8a_{4}v_{10}-fv_{6}^{2}}{8a_{4}v_{1}}\,,\qquad\tilde{M}^{(0)}_{11}\tilde{M}^{(0)}_{33}-\left(\tilde{M}^{(0)}_{13}\right)^{2}=0\,,$
(84)

under the use of relations among the coefficients given in Appendix A. We then
pick up the next-to leading order contribution proportional to $L^{0}$ in Eq.
(82). Although the resultant equation is complicated, we can simplify it by
using the following relations among coefficients appearing in Eq. (83),

$\displaystyle\frac{\tilde{K}^{(0)}_{23}}{\tilde{K}^{(0)}_{13}}=\frac{\tilde{K}^{(1)}_{23}}{\tilde{K}^{(1)}_{13}}=\frac{fb_{3}}{2ra_{4}}\,,\qquad\tilde{K}^{(0)}_{12}\tilde{M}^{(0)}_{33}-\tilde{K}^{(0)}_{23}\tilde{M}^{(0)}_{13}=rh\phi^{\prime}\tilde{K}_{22}^{(0)}\tilde{K}_{33}^{(0)}\,,$
$\displaystyle\tilde{K}^{(1)}_{11}\tilde{M}^{(0)}_{33}-2\tilde{K}^{(1)}_{13}\tilde{M}^{(0)}_{13}=r^{2}h^{2}\tilde{M}^{(0)}_{33}(\phi^{\prime
2}\tilde{K}_{22}^{(0)}+4r^{2}\alpha_{7})\,,$ (85)

together with Eq. (84). We also introduce the following two quantities,

$\displaystyle M_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{\tilde{M}_{11}^{(0)}{\tilde{M}}_{33}^{(1)}-2\tilde{M}_{13}^{(0)}{\tilde{M}}_{13}^{(1)}+\tilde{M}_{33}^{(0)}{\tilde{M}}_{11}^{(1)}}{4r^{2}fh^{2}\alpha_{7}\tilde{M}_{33}^{(0)}}\,,$
(86) $\displaystyle M_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{\phi^{\prime}(rh\phi^{\prime}{\tilde{M}}_{22}^{(0)}\tilde{M}_{33}^{(0)}+2\tilde{M}_{13}^{(0)}{\tilde{M}}_{23}^{(0)}-2{\tilde{M}}_{12}^{(0)}\tilde{M}_{33}^{(0)})}{4rfhv_{1}\alpha_{7}\tilde{M}_{33}^{(0)}}\,.$
(87)

Then, the next-to leading order contribution proportional to $L^{0}$ in Eq.
(82) can be expressed as

$\left(c_{\Omega}^{2}+\frac{r^{2}\tilde{M}_{33}^{(0)}}{f\tilde{K}_{33}^{(0)}}\right)\left[\left(c_{\Omega}^{2}+M_{1}+M_{2}\right)\left(c_{\Omega}^{2}+\frac{r^{2}\tilde{M}_{22}^{(0)}}{f\tilde{K}_{22}^{(0)}}\right)-\frac{r^{2}\left(4f\alpha_{7}M_{2}+\phi^{\prime
2}\tilde{M}_{22}^{(0)}\right)^{2}}{4f^{2}\phi^{\prime
2}\alpha_{7}\tilde{K}_{22}^{(0)}}\right]=0\,,$ (88)

where the matrix components appearing in the above expression and $M_{1}$,
$M_{2}$ are given by

$\displaystyle K_{22}^{(0)}=e_{1}\,,\qquad
M_{11}^{(0)}=-\frac{r^{2}v_{6}^{2}}{4v_{1}}\,,\qquad
M_{13}^{(0)}=-\frac{rv_{6}}{2}\,,\qquad M_{33}^{(0)}=-v_{1}\,,$ $\displaystyle
M_{11}^{(1)}=r^{2}d_{4}-\frac{m_{1}^{2}}{4a_{4}^{2}v_{9}}+\frac{rh(m_{1}-ra_{4}v_{8})m_{2}}{a_{4}v_{1}v_{6}}+\frac{1}{4}\left(\frac{rv_{6}m_{1}}{a_{4}v_{9}}-\frac{2r^{2}hm_{2}}{v_{1}}\right)^{\prime}\,,$
$\displaystyle
M_{33}^{(1)}=-\frac{m_{3}^{2}}{a_{4}^{2}v_{9}}-\frac{hA_{0}^{\prime}v_{1}m_{4}}{a_{4}}+\left(\frac{v_{1}m_{3}}{a_{4}v_{9}}\right)^{\prime}\,,$
$\displaystyle
M_{13}^{(1)}=-\frac{m_{1}m_{3}}{2a_{4}^{2}v_{9}}-\frac{rhA_{0}^{\prime}m_{2}}{2a_{4}}+\frac{h(m_{1}-ra_{4}v_{8})m_{4}}{2a_{4}v_{6}}+\frac{1}{4}\left(\frac{v_{1}m_{1}+rv_{6}m_{3}}{a_{4}v_{9}}-rhm_{4}\right)^{\prime}$
$\displaystyle
M_{22}^{(0)}=e_{4}+\frac{2hc_{4}a_{6}}{a_{4}}-\frac{m_{5}^{2}}{4a_{4}^{2}v_{9}}\,,$
$\displaystyle
M_{12}^{(0)}=-\frac{rd_{3}}{2}-\frac{hc_{4}(m_{1}-ra_{4}v_{8})}{a_{4}v_{6}}-\frac{rha_{6}m_{2}}{2a_{4}v_{1}}+\frac{rv_{5}v_{6}}{4v_{1}}+\frac{m_{1}m_{5}}{4a_{4}^{2}v_{9}}+\frac{1}{4}\left(2rhc_{4}+\frac{r(2d_{2}v_{1}-v_{4}v_{6})}{2v_{1}}-\frac{rv_{6}m_{5}}{2a_{4}v_{9}}\right)^{\prime}\,,$
$\displaystyle
M_{23}^{(0)}=\frac{hA_{0}^{\prime}c_{4}v_{1}}{a_{4}}+\frac{v_{5}}{2}+\frac{m_{3}m_{5}}{2a_{4}^{2}v_{9}}-\frac{ha_{6}m_{4}}{2a_{4}}-\frac{1}{4}\left(v_{4}+\frac{v_{1}m_{5}}{a_{4}v_{9}}\right)^{\prime}\,,$
(89)

with the shortcut notations

$\displaystyle m_{1}\equiv
ra_{4}v_{8}+\left(a_{4}+ra_{9}-\frac{rA_{0}^{\prime}v_{6}}{2}\right)v_{6}\,,\qquad
m_{2}\equiv
2c_{5}v_{1}+\left(A_{0}^{\prime}v_{1}+\frac{\phi^{\prime}v_{4}}{2}\right)v_{6}\,,$
$\displaystyle m_{3}\equiv
a_{4}v_{1}^{\prime}-\frac{A_{0}^{\prime}v_{1}v_{6}}{2}\,,\qquad m_{4}\equiv
2A_{0}^{\prime}v_{1}+\phi^{\prime}v_{4}\,,\qquad m_{5}\equiv
a_{6}v_{6}+a_{4}v_{13}\,.$ (90)

In a manner analogous to the radial propagation speeds, the angular
propagation speed of the vector DOF decouples from the other two propagation
speeds in Eq. (88). Denoting the vector propagating speed as $c_{\Omega
3{\rm,even}}$, it is given by

$c_{\Omega
3{\rm,even}}^{2}=-\frac{r^{2}\tilde{M}_{33}^{(0)}}{f\tilde{K}_{33}^{(0)}}=\frac{r^{2}\left(fv_{6}^{2}-8a_{4}v_{10}\right)}{8fa_{4}v_{1}}\geq
0\,,$ (91)

whose positivity is required for the Laplacian stability of the perturbation
$V$. This forms our third angular Laplacian stability condition, as indicated
in the subscript666This is referred as the third angular Laplacian stability
condition since it belongs to the vector DOF, which appears as the third one
in (64). For a more convenient way of narration, this one is discussed before
the first and second angular Laplacian stability conditions.. The remaining
two propagation speeds associated with the gravitational and scalar DOFs
generally couple with each other and can be expressed as

$\displaystyle
c_{\Omega\pm{\rm,even}}^{2}=\frac{1}{2}\left[-\left(M_{1}+M_{2}+\frac{r^{2}\tilde{M}_{22}^{(0)}}{f\tilde{K}_{22}^{(0)}}\right)\pm\sqrt{\left(M_{1}+M_{2}-\frac{r^{2}\tilde{M}_{22}^{(0)}}{f\tilde{K}_{22}^{(0)}}\right)^{2}+\frac{r^{2}\left(4f\alpha_{7}M_{2}+\phi^{\prime
2}\tilde{M}_{22}^{(0)}\right)^{2}}{f^{2}\phi^{\prime
2}\alpha_{7}\tilde{K}_{22}^{(0)}}}\right]$ (92)

For the absence of Laplacian instabilities associate with the perturbation
$\psi$ and $\delta\phi$, we require that the squared propagation speeds given
in Eq. (92) are real and non-negative. These conditions [with the assistance
of the Vieta’s formula and Eq.(88)], with the nearest two below referred as
the first and second angular Laplacian stability condition, respectively, are
translated to give

$\displaystyle\mathcal{M}_{1}\equiv-\left(M_{1}+M_{2}+\frac{r^{2}{\tilde{M}}_{22}^{(0)}}{f{\tilde{K}}_{22}^{(0)}}\right)\geq
0\,,$ (93)
$\displaystyle\mathcal{M}_{2}\equiv\frac{r^{2}(M_{1}+M_{2}){\tilde{M}}_{22}^{(0)}}{f{\tilde{K}}_{22}^{(0)}}-\frac{r^{2}\left(4f\alpha_{7}M_{2}+\phi^{\prime
2}{\tilde{M}}_{22}^{(0)}\right)^{2}}{4f^{2}\phi^{\prime
2}\alpha_{7}{\tilde{K}}_{22}^{(0)}}\geq 0\,,$ (94)

and

$\mathcal{M}_{1}^{2}-4\mathcal{M}_{2}=\left(M_{1}+M_{2}-\frac{r^{2}\tilde{M}_{22}^{(0)}}{f\tilde{K}_{22}^{(0)}}\right)^{2}+\frac{r^{2}\left(4f\alpha_{7}M_{2}+\phi^{\prime
2}\tilde{M}_{22}^{(0)}\right)^{2}}{f^{2}\phi^{\prime
2}\alpha_{7}\tilde{K}_{22}^{(0)}}\geq 0\,,$ (95)

with (95) guaranteeing that $c_{\Omega\pm{\rm,even}}^{2}$ are real. Imposing
the Laplacian stability condition (38) in the odd-parity sector, the quantity
$\alpha_{7}$ in the third condition can not be negative. Moreover, on using
the relations of coefficients in Appendix A and the definition of ${\cal
P}_{1}$ given in Eq. (74), the quantity $\tilde{K}_{22}^{(0)}$ in the third
condition can be written as

$\tilde{K}_{22}^{(0)}=e_{1}=\frac{\cal K}{\sqrt{fh}\phi^{\prime 2}}\,,$ (96)

which must be positive since ${\cal K}$ corresponds to the no-ghost condition
derived in Eq. (75). These facts show that the last condition (95) is
automatically guaranteed as long as the Laplacian stability in the odd-parity
sector and the no-ghost condition in the even-parity sector are satisfied.

In spite of the complexity of stability conditions (75), (81), (91), (93), and
(94), some of the stability conditions against even-parity perturbations can
be omitted since they are redundant with those against odd-parity
perturbations as we have already seen. Let us summarize that before closing
this subsection. For the absence of ghosts, we found that only one additional
condition, $\mathcal{K}>0$, to the odd mode stability is required where
$\mathcal{K}$ is given in Eq. (75). The tensor, scalar, and vector propagation
squared speeds along the radial direction, i.e., $c_{r1{\rm,even}}^{2}$,
$c_{r2{\rm,even}}^{2}$, and $c_{r3{\rm,even}}^{2}$, are given in Eqs. (79),
(80) and (78), respectively. Among them, we find that the tensor and vector
propagation speeds are equivalent to the radial propagation speeds of
corresponding DOFs against odd-parity perturbations given in Eqs. (34) and
(35), respectively. Hence, the radial Laplacian stability condition against
even-parity perturbations gives rise to only one additional condition
$c_{r2{\rm,even}}^{2}\geq 0$ to the odd mode stability. On the other hand, the
angular propagation speeds against even-parity perturbations are independent
of the odd mode stability. Hence, the angular Laplacian stabilities require
three conditions. The positivity of $c_{\Omega 3{\rm,even}}^{2}$ is satisfied
under the condition (91). The tensor and scalar propagation squared speeds,
$c_{\Omega\pm{\rm,even}}^{2}$, are guaranteed to be positive under the two
conditions (93) and (94).

In summary, we need to consider five additional conditions (75), (80), (91),
(93), and (94), to the stability of odd-parity sector. Nevertheless, it would
be still very hard to carry out the stability analysis for the most general
case. Instead, we shall process to next section and consider several concrete
models.

###  IV.3 Linear stability conditions for $l=0$

We consider the monopole perturbation characterized by $l=0$, i.e., $L=0$, in
which case the quantities $h_{0}$, $h_{1}$, and $G$ identically vanish away
from the second-order action of the even-parity sector Kase and Tsujikawa
(2023). This means that one can use the gauge DOFs on ${\cal T}$ and $\Theta$
for our purpose other than to eliminate $h_{0}$ and $G$ as in Eq. (45).
However, the gauge DOF will not be completely fixed in such a case. Thus, we
adopt the same gauge characterized by Eq. (45) as in the case for $l\geq 2$.
Substituting $L=0$ into the second-order action (49) with Eqs. (50) and (52),
it reduces to

$\displaystyle{\cal S}_{\rm even}^{l=0}$ $\displaystyle=$
$\displaystyle\int{\rm d}t\,{\rm d}r\Bigg{[}v_{1}\left\\{2V\left(\delta
A_{0}^{\prime}-\dot{\delta
A}_{1}+\frac{v_{3}H_{2}+v_{4}\delta\phi^{\prime}+v_{5}\delta\phi}{2v_{1}}\right)-V^{2}\right\\}-\frac{\left(v_{3}H_{2}+v_{4}\delta\phi^{\prime}+v_{5}\delta\phi\right)^{2}}{4v_{1}}$
(97)
$\displaystyle+\left(\Phi^{\prime}+A_{0}^{\prime}v_{1}V\right)H_{0}-\frac{2}{f}\dot{\Phi}H_{1}+\left(c_{2}\delta\phi^{\prime}+c_{3}\delta\phi\right)H_{2}+c_{6}H_{2}^{2}+e_{1}\dot{\delta\phi}^{2}+e_{2}\delta\phi^{\prime
2}+e_{3}\delta\phi^{2}\Bigg{]}\,,$

where we introduced

$\Phi\equiv
a_{1}\delta\phi^{\prime}+\left(a_{2}-a_{1}^{\prime}-\frac{1}{2}A_{0}^{\prime}v_{4}\right)\delta\phi+a_{3}H_{2}\,,$
(98)

and used the relations among coefficients in Appendix A. Compared to Eqs. (50)
and (52) for $l\neq 0$, the perturbations $H_{1}$, $\delta A_{0}$, $\delta
A_{1}$ do not possess the quadratic terms in Eq. (97). This fact shows that,
in addition to $H_{0}$, these variables also reduce to Lagrange multipliers
giving rise to constraints on the other variables for $l=0$. Indeed, the
variation of the action (97) with respect to $\delta A_{0}$, $\delta A_{1}$,
$H_{0}$, and $H_{1}$ lead to the following constraints

$\displaystyle(v_{1}V)^{\prime}=0\,,$ (99) $\displaystyle\dot{V}=0\,,$ (100)
$\displaystyle\Phi^{\prime}+A_{0}^{\prime}v_{1}V=0\,,$ (101)
$\displaystyle\dot{\Phi}=0\,,$ (102)

respectively. Integrating Eqs. (99) and (100), we obtain

$V=\frac{{\cal C}_{1}}{v_{1}}\,,$ (103)

where ${\cal C}_{1}$ denoting an integration constant. Substituting this
solution into Eq. (101) and integrating it together with Eq. (102), we obtain

$\Phi={\cal C}_{2}-{\cal C}_{1}A_{0}\,,$ (104)

where ${\cal C}_{2}$ is an integration constant. From the definition of $\Phi$
in Eq. (98), the above solution leads a constraint on the perturbation $H_{2}$
as

$H_{2}=\frac{1}{a_{3}}\left[{\cal C}_{2}-{\cal
C}_{1}A_{0}-a_{1}\delta\phi^{\prime}-\left(a_{2}-a_{1}^{\prime}-\frac{1}{2}A_{0}^{\prime}v_{4}\right)\delta\phi\right]\,.$
(105)

We substitute Eqs. (103), (104), and (105) into the second-order action (97)
for $l=0$. This operation removes all the metric perturbations, i.e., $H_{0}$,
$H_{1}$, and $H_{2}$, and the vector perturbation $V$ from the action. In
other words, we find that the monopole perturbation is governed by the scalar
perturbation $\delta\phi$ as a single DOF. Omitting the integration constants
${\cal C}_{1}$ and ${\cal C}_{2}$ irrelevant to the dynamics of perturbations,
the resultant action is expressed as

${\cal S}_{\rm even}^{l=0}=\int{\rm d}t\,{\rm
d}r\left[e_{1}\dot{\delta\phi}^{2}+\left(e_{2}-\frac{v_{4}^{2}}{4v_{1}}\right)\delta\phi^{\prime
2}+M_{0}\delta\phi^{2}\right]\,,$ (106)

where $M_{0}$ is given by

$\displaystyle M_{0}\equiv
e_{3}+\frac{b_{3}(b_{3}c_{6}-b_{4}c_{3})}{b_{4}^{2}}-\frac{(b_{3}v_{3}-b_{4}v_{5})^{2}}{4b_{4}^{2}v_{1}}+\frac{1}{2}\left(\frac{b_{3}c_{2}}{b_{4}}-\frac{(b_{3}v_{3}-b_{4}v_{5})v_{4}}{2b_{4}v_{1}}\right)^{\prime}\,.$
(107)

The above action shows that the no-ghost condition for the monopole
perturbation is given by $e_{1}>0$. As we have seen in Eq. (96), this is
equivalent to the no-ghost condition for the perturbations with $l\geq 2$. The
propagation speed square along the radial direction, $c_{r,{\rm
even}}^{2}=-(e_{2}-v_{4}^{2}/v_{1})/e_{1}$, also coincides with that for
$l\geq 2$ derived in Eq. (80). Thus, we have shown that the stability
conditions derived for $l\geq 2$ also ensure the absence of instabilities in
the monopole perturbation.

###  IV.4 Linear stability conditions for $l=1$

We proceed to study the dipole perturbations characterized by $l=1$ ($L=2$) in
which case the metric perturbations $K$ and $G$ appear in the second-order
action only through the combination of the form $G-K$ Kase and Tsujikawa
(2023). This means that, instead of using the gauge DOFs of ${\cal R}$ and
$\Theta$ to eliminate both $G$ and $K$ as we adopted in Eq. (45), it is
possible to keep one of them by eliminating the combination $G-K$ and fixing
either ${\cal R}$ or $\Theta$. Since Eq. (44) shows that the complete gauge
choice of $\Theta$ can be obtained only via the transformation of the
perturbation $G$, we use this gauge DOF to eliminate $G$. The residual gauge
DOF ${\cal R}$ can also be completely fixed via the transformation of
$\delta\phi$. Thus, we choose the following gauge for the dipole perturbation,

$h_{0}=0\,,\qquad\delta\phi=0\,,\qquad G=K\,.$ (108)

The elimination of non-dynamical variables in the second-order action for
$l=1$ can be operated in the same way for $l\geq 2$ which we discussed below
Eq. (62). After this process, the resultant action is composed of two
dynamical variables, $\psi$ and $V$, showing that the dipole perturbation
possesses one less propagating DOF than the case with $l\geq 2$. In an
analogous way to Sec. IV.2.1, we find that the no-ghost conditions for these
two DOFs coincide with Eqs. (70) and (72) substituted $L=2$. The squared
propagation speeds along the radial direction for $l=1$ coincide with Eqs.
(78) and (80), i.e., the propagation speeds of vector and scalar field
perturbations for $l\geq 2$, respectively. These facts show that the dipole
perturbation is governed by vector and scalar field perturbations.
Consequently, we find that the stability of the dipole perturbation does not
give rise to additional conditions to those derived for $l\geq 2$.

##  V Application to the concrete models

In Ref. Heisenberg and Tsujikawa (2018), three different concrete models
possessing hairy BH solutions on the spherically symmetric spacetime were
proposed based on the $U(1)$ GI SVT theories. The stability of such models
against odd-parity perturbations are studied in Ref. Heisenberg  _et al._
(2018). These models are characterized by the following choices of functions

$\displaystyle{\rm Model\ 1}:\quad$ $\displaystyle f_{3}=\beta_{3},\qquad
f_{4}=0\,,$ (109) $\displaystyle{\rm Model\ 2}:\quad$ $\displaystyle
f_{3}=\beta_{3},\qquad f_{4}=\beta_{4}\,,$ (110) $\displaystyle{\rm Model\
3}:\quad$ $\displaystyle f_{3}=\beta_{3},\qquad f_{4}=\beta_{4}X\,,$ (111)

respectively, where $\beta_{3}$ and $\beta_{4}$ are arbitrary constants under
the constraints led by the odd-parity stability analysis (and later will be
further confined by the even-parity ones). The other functions are common to
all three models and are chosen as $f_{2}=X+F$ plus
${\tilde{f}}_{3}={\tilde{f}}_{4}=0$. In the following, we shall discuss the
stability conditions against odd-parity perturbations, Eqs. (33) and (38), and
those against even-parity perturbations, Eqs. (75), (81), (91), (93) and (94),
for each model. Notice that, for the convenience of readers, some of the main
results are summarized at the end of this section in Table 1. One can move to
there directly when the concluding remarks are demanded.

###  V.1 Stability analysis for the Model 1

In this model, the quantities $\alpha_{1}$, $\alpha_{6}$, and $\alpha_{7}$
associated with stability conditions against odd-parity perturbations reduce
to

$\alpha_{1}=\frac{M_{\rm
pl}^{2}\sqrt{h}}{4\sqrt{f}}\,,\qquad\alpha_{6}=-\frac{\sqrt{fh}M_{\rm
pl}^{2}}{4r^{2}}\,,\qquad\alpha_{7}=\frac{M_{\rm pl}^{2}}{4r^{2}\sqrt{fh}}\,.$
(112)

This shows that one of the no-ghost conditions, $\alpha_{6}<0$, is trivially
satisfied in Eq. (33). Moreover, the propagation speed of gravitational DOF
along the radial and angular direction reduce to that of light, i.e.,
$c_{r1,{\rm odd}}^{2}=1=c_{\Omega 1,{\rm odd}}^{2}$. This fact shows that the
Laplacian stability conditions for the gravitational DOF in the odd-parity
sector are also automatically satisfied. The stability conditions associated
with the vector field perturbation remain nontrivial. Thus, in Eqs. (33) and
(38), the nontrivial stability conditions for this model are

$\alpha_{4}>0\,,\qquad\alpha_{2}^{2}-\alpha_{1}\alpha_{5}\geq
0\,,\qquad\alpha_{8}^{2}-4\alpha_{7}\alpha_{9}\geq 0\,,$ (113)

which are consistent with the Eqs. (3.23), (3.25) and (3.28) in Ref.
Heisenberg  _et al._ (2018). The above constraints can be numerically
translated to the phase space of $\\{\beta_{3},\beta_{4}\\}$ in the models
(110) and (111). An overall constraint on $\\{\beta_{3},\beta_{4}\\}$ from the
odd-parity sector can be found in Ref. Heisenberg  _et al._ (2018) by
comprehensively considering all the three models discussed in here. We shall
go back to this point in the following subsections.

In the eve-parity sector, the no-ghost condition represented by Eq. (75)
reduces to

$\displaystyle\mathcal{K}$ $\displaystyle=$
$\displaystyle\frac{2r_{h}^{4}h^{2}A_{0}^{\prime
4}{\tilde{\beta}}_{3}^{2}}{M_{\rm pl}^{2}f^{2}}\,,$ (114)

where we have defined the dimensionless factor ${\tilde{\beta}}_{3}\equiv
M_{\rm pl}\beta_{3}/r_{h}^{2}$ and $r_{h}$ stands for the radii of the event
horizon so that we have $f(r_{h})=h(r_{h})=0$. Since Eq. (114) shows that the
quantity $\mathcal{K}$ is positive for any non-zero value of
${\tilde{\beta}}_{3}$, the third no-ghost condition will always get
satisfied777 According to the (3.14) and (3.15) of Heisenberg and Tsujikawa
(2018), we will have $f/h\to 1$ as $r\to r_{h}$. Thus, $\mathcal{K}$ won’t be
zero at the event horizon so that it has to stay positive. In addition, notice
that, here we are assuming $A_{0}^{\prime}\neq 0$.. We note that the
background equations (12)-(15) have been used during simplifications and we
shall apply the same operation without additional alerts in the following.
Regarding Eqs. (14)-(15), we especially notice that the conditions
$J_{\phi}(r)=0$ and $J_{A}(r)=\text{constant}\equiv r_{h}^{2}\kappa$ hold as
discussed in the end of Sec. II. Thus, we have

$\displaystyle\phi^{\prime}=\frac{2h\tilde{\beta}_{3}A_{0}^{\prime
2}r_{h}^{2}}{rfM_{\text{pl}}}\,,\qquad
A_{0}^{\prime}=\frac{\kappa\sqrt{f}r_{h}^{2}M_{\text{pl}}}{r\sqrt{h}\left(4h\tilde{\beta}_{3}r_{h}^{2}\phi^{\prime}+rM_{\text{pl}}\right)}\,.$
(115)

Here, $\kappa$ is temporarily borrowed to denote a dimensionless constant.
With the polynomial solutions given by Eqs. (3.14)-(3.16) in Ref. Heisenberg
and Tsujikawa (2018) and evaluation of (18) at $r_{h}$, we notice that for
Model 1 we have $\kappa=\sqrt{2\mu}M_{\text{pl}}/r_{h}$, where $\mu\in(0,1)$
so that $\left(r_{h}/M_{\text{pl}}\right)\kappa\in(0,\sqrt{2})$ Heisenberg
_et al._ (2018).

On the other hand, for Model 1, the squared second propagation speed given by
Eq. (80) reduces to

$\displaystyle c_{r2{\rm,even}}^{2}=1\,,$ (116)

which shows that the second Laplacian stability condition along the radial
direction is automatically satisfied.

Regarding the third Laplacian stability condition along the angular direction,
we plug the full expressions of $v_{6}$, $v_{10}$, etc., into Eq. (91), we
find that

$\displaystyle c_{\Omega 3{\rm,even}}^{2}$ $\displaystyle=$
$\displaystyle{\tilde{K}}_{33}\frac{2hL\left[2f^{2}h(h+1)\tilde{\beta}_{3}^{2}A_{0}^{\prime
2}r_{h}^{4}M_{\text{pl}}^{4}-fh^{2}r^{2}\tilde{\beta}_{3}^{2}A_{0}^{\prime
4}r_{h}^{4}M_{\text{pl}}^{2}-12h^{4}\tilde{\beta}_{3}^{4}A_{0}^{\prime
6}r_{h}^{8}+f^{3}r^{2}M_{\text{pl}}^{6}\right]^{2}}{r^{2}(fh)^{3/2}M_{\text{pl}}^{6}\left(8h^{2}\tilde{\beta}_{3}^{2}A_{0}^{\prime
2}r_{h}^{4}+fr^{2}M_{\text{pl}}^{2}\right)^{3}}\,.$ (117)

From its mathematical structure we notice that, $c_{\Omega 3{\rm,even}}^{2}$
is definitely a non-negative quantity provided that $\tilde{K}_{33}>0$ holds
under the stability of odd-parity sector. Thus, this angular Laplacian
stability condition gets satisfied.

We proceed to the first angular Laplacian stability condition in which the
story becomes a little sophisticated. To deal with this problem, we first
convert $\mathcal{M}_{1}$ given in Eq. (93) into a dimensionless form by (i)
solving Eqs. (115) for $A_{0}^{\prime}$ 888Notice that, this set of algebraic
equations have 3 roots originally. However, only one of them is real so we
keep it as our solution for $A_{0}^{\prime}$. and substituting it into
$\mathcal{M}_{1}$; (ii) replacing the model parameter $\beta_{3}$ in
$\mathcal{M}_{1}$ with the dimensionless ${\tilde{\beta}}_{3}$; (iii) changing
the variable $r$ to the dimensionless one $\xi\equiv r_{h}/r$; (iv) taking
advantage of the freedom for choosing unit system, setting
$M_{\text{pl}}=r_{h}$ (cf., footnote #1). As a result, a dimensionless form of
$\mathcal{M}_{1}$ is obtained as a function of
$\\{f,h,\xi,\kappa,{\tilde{\beta}}_{3}\\}$. We postpone its lengthy explicit
expression to Appendix B. Knowing the fact that $f,h\in[0,1]$, $\xi\in(0,1]$
since we have $r\in[r_{h},+\infty)$, $\kappa\in(0,\sqrt{2})$ and
$|{\tilde{\beta}}_{3}|\lesssim{\cal O}(1)$Heisenberg  _et al._ (2018), we can
attempt to determine the range of $\mathcal{M}_{1}$ in the 5-dimension phase
space spanned by $\\{f,h,\xi,\kappa,{\tilde{\beta}}_{3}\\}$. In practice, all
of these could be fulfilled numerically on Mathematica. It turns out that the
minimum of $\mathcal{M}_{1}$ in this phase space is indeed positive, which is
about999Strictly speaking, the true minimum should be bigger than this, since
$f$ and $h$ are correlated so that certain regions of this phase space won’t
be reached. $2.001$. This result guarantees this angular Laplacian stability
condition.

The second angular Laplacian stability condition (94) is even more
complicated. To analyze $\mathcal{M}_{2}$, new strategies have to be selected.
After some attempts, we know a posteriori that new constraints on
${\tilde{\beta}}_{3}$ have to be performed for this stability condition to
hold. Now the mission is to determine the new constraint. For this goal and by
referring the analysis to the first angular Laplacian stability condition in
which the minimum of $\mathcal{M}_{1}$ is achieved near the $r=r_{h}$, as well
as the stability conditions in odd-parity sector studied in Ref. Heisenberg
_et al._ (2018) in which the previously mentioned constraint on
${\tilde{\beta}}_{3}$ is obtained near the $r=r_{h}$, we shall focus on the
behavior of $\mathcal{M}_{2}$ near the horizon $r\to r_{h}$. As usual, a
dimensionless form of $\mathcal{M}_{2}$ could be obtained by following the
same procedure mentioned earlier. However, since we are focusing on the $r\to
r_{h}$ limit, the phase space now reduces to 2 dimension, spanned by
$\\{\mu,{\tilde{\beta}}_{3}\\}$. Now the $\mathcal{M}_{2}$ could be written as

$\displaystyle\left.\mathcal{M}_{2}\right|_{r\to r_{h}}$ $\displaystyle=$
$\displaystyle\Big{[}6144(\mu-1)^{2}\mu^{5}(3\mu-4)\tilde{\beta}_{3}^{8}-3200(\mu-1)^{2}\mu^{4}\tilde{\beta}_{3}^{6}+16\mu^{2}\left(9\mu^{2}-10\mu+1\right)\tilde{\beta}_{3}^{4}$
(118)
$\displaystyle-8192\mu^{6}\left(2\mu^{2}-5\mu+3\right)^{2}\tilde{\beta}_{3}^{10}-8(\mu-1)\mu\tilde{\beta}_{3}^{2}+1\Big{]}{\left[1-4(\mu-1)\mu\tilde{\beta}_{3}^{2}\right]^{-2}}\,.$

We plot it out in Fig. 1 of Appendix C. By looking at the contour of the
allowed region for $\mathcal{M}_{2}$, we can read off the new upper bound for
legal $|{\tilde{\beta}}_{3}|$, which is approximately given by
$|{\tilde{\beta}}_{3}|\lesssim 0.441$. In addition, it’s worth mentioning here
that, by adopting this new upper bound and the numerical method for analyzing
$\mathcal{M}_{1}$ with Mathematica, the minimum of $\mathcal{M}_{2}$ in the
allowed region of phase space $\\{f,h,\xi,\kappa,{\tilde{\beta}}_{3}\\}$ is
indeed found to be positive. Thus, we conclude that, under the new constraint
$|{\tilde{\beta}}_{3}|\lesssim 0.441$, the second angular Laplacian stability
condition is guaranteed.

###  V.2 Stability analysis for the Model 2

As before, we first consider the third no-ghost condition (75). For Model 2
characterized by Eq. (110), $\mathcal{K}$ becomes

$\displaystyle\mathcal{K}$ $\displaystyle=$
$\displaystyle\frac{2h^{2}\tilde{\beta}_{3}^{2}A_{0}^{\prime
4}r_{h}^{4}}{f^{2}M_{\text{pl}}^{2}}\,,$ (119)

which happens to be identical to Eq. (114). Thus, for any non-zero
${\tilde{\beta}}_{3}$, the third no-ghost condition will always get satisfied.
As mentioned in last subsection, the background equations (12)-(15) have been
used during simplifications. By the same procedure to obtain Eq. (115), Eqs.
(14)-(15) lead to

$\displaystyle\phi^{\prime}=\frac{2h\tilde{\beta}_{3}A_{0}^{\prime
2}r_{h}^{2}}{frM_{\text{pl}}}\,,\qquad
A_{0}^{\prime}=\frac{\kappa\sqrt{f}r_{h}^{2}M_{\text{pl}}}{\sqrt{h}\left(-8h\tilde{\beta}_{4}r_{h}^{2}M_{\text{pl}}+8\tilde{\beta}_{4}r_{h}^{2}M_{\text{pl}}+4hr\tilde{\beta}_{3}r_{h}^{2}\phi^{\prime}+r^{2}M_{\text{pl}}\right)}\,,$
(120)

where we have defined the dimensionless factor ${\tilde{\beta}}_{4}\equiv
r_{h}^{-2}\beta_{4}$. With the polynomial solutions given in Eqs. (4.5)-(4.7)
of Ref. Heisenberg and Tsujikawa (2018) and evaluate Eq. (18) at $r_{h}$, we
notice that for Model 2 we have
$\kappa=\sqrt{2\mu(1+{\tilde{\beta}}_{4})}M_{\text{pl}}/r_{h}$, where
$\mu\in(0,1)$. Since the odd-mode stability conditions give the constraint
${\tilde{\beta}}_{4}\in(-0.125,0.251)$ (see also Eqs. (4.12)-(4.14), (4.20)
and (4.21) in Ref. Heisenberg  _et al._ (2018)), we obtain that
$\left(r_{h}/M_{\text{pl}}\right)\kappa\in(0,1.59)$. To manifest the
constraints on ${\tilde{\beta}}_{4}$, we plot out the allowed region given by
the stability conditions against odd-parity perturbations, i.e., Eqs.
(4.12)-(4.14), (4.20) and (4.21) in Ref. Heisenberg  _et al._ (2018) in the
3-dimension phase space $\\{{\tilde{\beta}}_{3},{\tilde{\beta}}_{4},\mu\\}$ in
the panel (a) of Fig. 2 in Appendix C, where the constraints $\mu\in(0,1)$ and
$|{\tilde{\beta}}_{3}|\lesssim 0.441$ are assumed. From this figure we learn
that, to meet these constrains for any legal arbitrary $\mu$, there must be
upper and lower bounds on ${\tilde{\beta}}_{4}$. We further observe that, the
most stringent constraints are from the $\mu\approx 0$ case. Thus, we plot out
the allowed region in $\\{{\tilde{\beta}}_{3},{\tilde{\beta}}_{4}\\}$ phase
space in the panel (b) of Fig. 2 in Appendix C by setting $\mu=0$. Finally,
from this panel we read off the constraints of ${\tilde{\beta}}_{4}$, which is
about ${\tilde{\beta}}_{4}\in(-0.125,0.251)$101010Interestingly, an upper
bound to ${\tilde{\beta}}_{4}$, which is quite similar to the one found in
here, could also be obtained by solely doing the even-parity stability
analysis. Ignoring the tolerable numerical errors, we conclude that the
results from odd- and even-parity sectors are consistent. .

On the other hand, for Model 2, the squared second radial propagation speed
given in Eq. (80) reduces to

$\displaystyle c_{r2{\rm,even}}^{2}=1\,,$ (121)

which shows that the second radial Laplacian stability condition is
automatically satisfied.

Next, let us substitute the full expressions of $v_{6}$, $v_{10}$, etc., into
Eq. (91) to obtain

$\displaystyle c_{\Omega 3{\rm,even}}^{2}$ $\displaystyle=$
$\displaystyle\frac{{\tilde{K}}_{33}}{r^{2}(fh)^{3/2}M_{\text{pl}}^{6}\left[8h^{2}\tilde{\beta}_{3}^{2}A_{0}^{\prime
2}r_{h}^{4}+fM_{\text{pl}}^{2}\left(8\tilde{\beta}_{4}r_{h}^{2}(1-h)+r^{2}\right)\right]^{3}}\,$
$\displaystyle\times 2hL\Big{\\{}2f^{2}hA_{0}^{\prime
2}r_{h}^{2}M_{\text{pl}}^{4}\left[\left(\tilde{\beta}_{3}^{2}+8\tilde{\beta}_{4}^{2}\right)r_{h}^{2}+h\left(\tilde{\beta}_{3}^{2}-8\tilde{\beta}_{4}^{2}\right)r_{h}^{2}+r^{2}\tilde{\beta}_{4}\right]$
$\displaystyle-fh^{2}\tilde{\beta}_{3}^{2}A_{0}^{\prime
4}r_{h}^{4}M_{\text{pl}}^{2}\left(-32h\tilde{\beta}_{4}r_{h}^{2}+8\tilde{\beta}_{4}r_{h}^{2}+r^{2}\right)-12h^{4}\tilde{\beta}_{3}^{4}A_{0}^{\prime
6}r_{h}^{8}+f^{3}M_{\text{pl}}^{6}\left(4(h-1)\tilde{\beta}_{4}r_{h}^{2}+r^{2}\right)\Big{\\}}^{2}.$

Similar to Eq. (117), this result shows that $c_{\Omega 3{\rm,even}}^{2}$ is
definitely non-negative as long as $\tilde{K}_{33}>0$ holds, so that this
angular Laplacian stability is guaranteed.

Moving to the first angular Laplacian stability condition, the story is even
more sophisticated than that of Model 1 due to the presence of
${\tilde{\beta}}_{4}$. Because of this, the analysis of $\mathcal{M}_{1}$ will
be divided into two parts. The first attempt is about the behavior around
$r_{h}$. As usual, a dimensionless form of $\mathcal{M}_{1}$ is needed. The
basic steps for obtaining that is the same as what we did in last subsection.
The only difference is that the resultant $\mathcal{M}_{1}$ is now a function
of six variables, i.e.,
$\\{f,h,\xi,\kappa,{\tilde{\beta}}_{3},{\tilde{\beta}}_{4}\\}$. On top of
that, it became more accessible to expand $\mathcal{M}_{1}$ around the metric
horizon. To the lowest order of $(r-r_{h})$, that leads to

$\displaystyle\left.\mathcal{M}_{1}\right|_{r\to r_{h}}$ $\displaystyle=$
$\displaystyle
2\Big{\\{}4\kappa^{10}\tilde{\beta}_{3}^{4}\tilde{\beta}_{4}+\kappa^{8}\tilde{\beta}_{3}^{4}\left(-128\tilde{\beta}_{4}^{2}+8\tilde{\beta}_{4}+3\right)$
(123)
$\displaystyle-2\kappa^{6}\left(8\tilde{\beta}_{4}+1\right)^{2}\left[\left(3-8\tilde{\beta}_{4}\right)\tilde{\beta}_{3}^{4}+2\tilde{\beta}_{4}\left(64\tilde{\beta}_{4}^{2}+48\tilde{\beta}_{4}+5\right)\tilde{\beta}_{3}^{2}+8\tilde{\beta}_{4}^{3}\left(8\tilde{\beta}_{4}+1\right)^{2}\right]$
$\displaystyle-2\kappa^{2}\left(8\tilde{\beta}_{4}+1\right)^{6}\left[\left(8\tilde{\beta}_{4}+1\right)\tilde{\beta}_{3}^{2}+2\tilde{\beta}_{4}\left(80\tilde{\beta}_{4}^{2}+14\tilde{\beta}_{4}-1\right)\right]$
$\displaystyle+4\tilde{\beta}_{4}\left[2\left(16\tilde{\beta}_{4}+7\right)\tilde{\beta}_{3}^{2}+\tilde{\beta}_{4}\left(256\tilde{\beta}_{4}^{2}+88\tilde{\beta}_{4}+7\right)\right]\left(8\kappa\tilde{\beta}_{4}+\kappa\right)^{4}+\left(4\tilde{\beta}_{4}-1\right)\left(8\tilde{\beta}_{4}+1\right)^{9}\Big{\\}}$
$\displaystyle\times\left(8\tilde{\beta}_{4}+1\right)^{-4}\left(-4\left(\kappa^{2}-4\right)\tilde{\beta}_{4}+64\tilde{\beta}_{4}^{2}+1\right)^{-1}$
$\displaystyle\times\left[\kappa^{4}\tilde{\beta}_{3}^{2}-2\kappa^{2}\left(8\tilde{\beta}_{4}+1\right)\left(\tilde{\beta}_{3}^{2}+8\tilde{\beta}_{4}^{2}+\tilde{\beta}_{4}\right)+\left(4\tilde{\beta}_{4}-1\right)\left(8\tilde{\beta}_{4}+1\right)^{3}\right]^{-1}\,.$

On using this expression, we can determine the minimum of $\mathcal{M}_{1}$ in
the phase space of $\\{\mu,{\tilde{\beta}_{3}},{\tilde{\beta}_{4}}\\}$ with
the built-in functions of Mathematica. Keeping in mind that $\mu\in(0,1)$,
${\tilde{\beta}_{3}}\in(-0.441,0.441)$ and
${\tilde{\beta}_{4}}\in(-0.125,0.251)$, it turns out that the minimum of
$\mathcal{M}_{1}$ is a positive number, which supports this angular Laplacian
stability condition for Model 2. A similar treatment to $\mathcal{M}_{2}$
brings us a result in the form of Eq. (123). Because of its tedious
expression, we put that in Appendix B. The result shows that the minimum of
$\mathcal{M}_{2}$ in the phase space of
$\\{\mu,{\tilde{\beta}_{3}},{\tilde{\beta}_{4}}\\}$ is also determined to be a
positive number. Therefore, it seems like the currently known constraints are
sufficient to support this angular Laplacian stability condition.

We note that the above discussion is based on just an analysis at the $r\to
r_{h}$ limit. That brings us to the second part of analysis. To determine the
ranges of $\mathcal{M}_{1,2}$ in the whole exterior space
$r\in[r_{h},+\infty)$, we shall mimic Heisenberg  _et al._ (2018) and plot
them out for chosen parameters. To do so, one has to first solve for the
background quantities $f$ and $h$ on the exterior space with the algebraic
expressions (120) of $\phi^{\prime}$ and $A_{0}^{\prime}$ in terms of $f$, $h$
as well as their derivatives. They are achieved by numerically integrating
Eqs. (12) and (13) from $\xi=1-\epsilon_{1}$ to $\xi=\epsilon_{2}$, where
factors $\epsilon_{1}$ and $\epsilon_{2}$ are chosen to be small enough to
cover the desired region (see, e.g., Ref. Zhang  _et al._ (2020) for more
details about the relative techniques used in doing the numerical integration
and searching for background solutions). Of course, in practice we were using
the dimensionless forms of Eqs. (12) and (13), and a new quantity similar to
the $\mathfrak{Z}(\xi)$ defined in Eq. (133) in Appendix B is introduced to
simplify the expressions. Since this kind of treatment is straightforward and
we have already briefly described the steps for similar scenarios in last
subsection, we omit the details here. We note that the Taylor expansion

$\displaystyle{\mathfrak{f}}(\xi=1-\epsilon_{1})=\left.{\mathfrak{f}}\right|_{\xi=1}+\frac{1}{1!}\left.\frac{d{\mathfrak{f}}}{d\xi}\right|_{\xi=1}(-\epsilon_{1})+\frac{1}{2!}\left.\frac{d^{2}{\mathfrak{f}}}{d\xi^{2}}\right|_{\xi=1}(-\epsilon_{1})^{2}\,,$
(124)

where ${\mathfrak{f}}$ stands for $f$ and $h$, was applied in generating the
boundary conditions at $\xi=1-\epsilon_{1}$ keeping in mind that
${\mathfrak{f}}(\xi=1)=0$. The derivatives appearing in Eq. (124) are obtained
by using Eqs. (12), (13) as well as their derivatives. It’s worth mentioning
here that, the polynomial solutions given in Ref. Heisenberg and Tsujikawa
(2018) were not utilized directly in generating the boundary conditions since
there is a remaining parameter corresponding to an overall factor of $f$,
which is unknown before obtaining the numerical solutions. This factor will be
fixed by normalizing $f$ according to the requirement $f\to 1$ at the spatial
infinity.

The solution for $f$ and $h$ are shown in the panel (a) of Fig. 3 in Appendix
C. Since the upper limit of ${\tilde{\beta}_{3}}$ is updated to $\approx
0.441$ compared to Ref. Heisenberg  _et al._ (2018), we adopted it with
${\tilde{\beta}_{4}}=0.1$ and $\mu=0.5$. We can find deviation between Model 2
and that of GR. With the solutions of $f$ and $h$ in hand, together with the
previously mentioned algebraic expressions for $A_{0}^{\prime}$ and
$\phi^{\prime}$, all the background quantities are known. Inserting them into
${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$, we further obtain their numerical
values. This part of results are exhibited in the panel (b) of Fig. 3 in which
${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ are plotted as functions of
$r_{h}/r$ instead of $r/r_{h}$ to show a larger scope. Fig. 3 shows that, both
these two quantities stay positive in the whole space outside the horizon.
This result further supports that the first and second angular Laplacian
stability conditions for Model 2111111 An ideal treatment is to test the
positiveness of ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ using the Monte
Carlo method and try to cover the whole allowed region in
$\\{\mu,{\tilde{\beta}_{3}},{\tilde{\beta}_{4}},\xi\\}$ phase space.
Nevertheless, as seen from, e.g., Eqs. (132)-(134), the complexity of our
problem has already made the calculations quite time-consuming and it is
preventing us from running a more complete analysis. In principle, much more
computational resources are needed for executing such an analysis. .

###  V.3 Stability analysis for the Model 3

We first consider the third no-ghost condition (75). For Model 3,
$\mathcal{K}$ reduces to

$\displaystyle\mathcal{K}$ $\displaystyle=$
$\displaystyle\frac{2hr^{2}\tilde{\beta}_{3}^{2}A_{0}^{\prime
4}r_{h}^{4}\left[4(1-h)\tilde{\beta}_{4}A_{0}^{\prime
2}\sqrt{fh^{7}}r_{h}^{2}+f^{3/2}h^{5/2}r^{2}\right]}{(fh)^{3/2}M_{\text{pl}}^{2}\left(6h^{2}\tilde{\beta}_{4}A_{0}^{\prime
2}r_{h}^{2}-4h\tilde{\beta}_{4}A_{0}^{\prime
2}r_{h}^{2}-fr^{2}\right)^{2}}\,,$ (125)

Thus, by looking at both the numerator and denominator of the above
expression, for any non-zero ${\tilde{\beta}}_{3}$, the third no-ghost
condition will always get satisfied. As mentioned in previous subsections, the
background equations (14)-(15) lead to

$\displaystyle\phi^{\prime}$ $\displaystyle=$
$\displaystyle\frac{2hr\tilde{\beta}_{3}A_{0}^{\prime
2}r_{h}^{2}}{M_{\text{pl}}\left(-6h^{2}\bar{\beta}_{4}A_{0}^{\prime
2}r_{h}^{2}+4h\bar{\beta}_{4}A_{0}^{\prime 2}r_{h}^{2}+fr^{2}\right)}\,,$
$\displaystyle A_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\kappa\sqrt{f}r_{h}^{2}M_{\text{pl}}}{\sqrt{h}\left(6h^{2}\tilde{\beta}_{4}r_{h}^{2}M_{\text{pl}}\phi^{\prime
2}-4h\tilde{\beta}_{4}r_{h}^{2}M_{\text{pl}}\phi^{\prime
2}+4hr\tilde{\beta}_{3}r_{h}^{2}\phi^{\prime}+r^{2}M_{\text{pl}}\right)}\,.$
(126)

With the polynomial solutions given by Eqs. (4.15)-(4.17) in Ref. Heisenberg
and Tsujikawa (2018) and evaluation of (18) at $r_{h}$, we notice that for
Model 3 we have $\kappa=\sqrt{2\mu}M_{\text{pl}}/r_{h}$, where $\mu\in(0,1)$.
It is quite necessary for one to notice that, the system given by Eq. (126)
forms a fifth-order algebraic equation of $A_{0}^{\prime}$. As a result, there
is no analytic expression of $\phi^{\prime}$ or $A_{0}^{\prime}$ in terms of
$f$, $h$, etc. according to the well-known Abel-Ruffini theorem. This is one
of the key factors which makes the Model 3 much intricate than that of Model
2.

On the other hand, for Model 3, the squared second radial propagation speed
given by Eq. (80) becomes

$\displaystyle c_{r2{\rm,even}}^{2}$ $\displaystyle=$
$\displaystyle\frac{3h^{2}\tilde{\beta}_{4}A_{0}^{\prime
2}r_{h}^{2}\phi^{\prime
2}+fM_{\text{pl}}^{2}}{2h^{2}\tilde{\beta}_{4}A_{0}^{\prime
2}r_{h}^{2}\phi^{\prime 2}+fM_{\text{pl}}^{2}}\,,$ (127)

which shows that the second radial Laplacian stability condition holds.

Next, the angular propagation speed given in Eq. (91) for the Model 3 is

$\displaystyle c_{\Omega 3{\rm,even}}^{2}$ $\displaystyle=$
$\displaystyle\frac{{{\tilde{K}}_{33}}}{M_{\text{pl}}^{6}r^{4}\sqrt{fh}\left(2h(2-3h)\bar{\beta}_{4}A_{0}^{\prime
2}r_{h}^{2}+fr^{2}\right)^{2}}$ (128) $\displaystyle\times
2L\Big{\\{}16h^{5}(3h-2)r^{2}r_{h}^{10}\bar{\beta}_{4}^{2}\tilde{\beta}_{3}^{2}\left(2\bar{\beta}_{4}M_{\text{pl}}^{2}+h\left(2r_{h}^{2}\tilde{\beta}_{3}^{2}-3M_{\text{pl}}^{2}\bar{\beta}_{4}\right)\right)A_{0}^{\prime
10}$
$\displaystyle+4fh^{4}r_{h}^{8}\bar{\beta}_{4}\Big{[}108h^{4}\bar{\beta}_{4}^{2}\left(3M_{\text{pl}}^{2}\bar{\beta}_{4}-2r_{h}^{2}\tilde{\beta}_{3}^{2}\right)M_{\text{pl}}^{4}-72h^{3}\bar{\beta}_{4}^{2}\left(12M_{\text{pl}}^{2}\bar{\beta}_{4}-5r_{h}^{2}\tilde{\beta}_{3}^{2}\right)M_{\text{pl}}^{4}$
$\displaystyle+4\bar{\beta}_{4}\left(8M_{\text{pl}}^{2}\bar{\beta}_{4}\left(2\bar{\beta}_{4}M_{\text{pl}}^{2}+r_{h}^{2}\tilde{\beta}_{3}^{2}\right)-3r^{4}\tilde{\beta}_{3}^{2}\right)M_{\text{pl}}^{2}$
$\displaystyle-8h\left(\left(2r_{h}^{2}\tilde{\beta}_{3}^{4}-3M_{\text{pl}}^{2}\bar{\beta}_{4}\tilde{\beta}_{3}^{2}\right)r^{4}+8M_{\text{pl}}^{4}\bar{\beta}_{4}^{2}\left(6\bar{\beta}_{4}M_{\text{pl}}^{2}+r_{h}^{2}\tilde{\beta}_{3}^{2}\right)\right)$
$\displaystyle+3h^{2}\left(\left(2r_{h}^{2}\tilde{\beta}_{3}^{4}-3M_{\text{pl}}^{2}\bar{\beta}_{4}\tilde{\beta}_{3}^{2}\right)r^{4}+8M_{\text{pl}}^{4}\bar{\beta}_{4}^{2}\left(36M_{\text{pl}}^{2}\bar{\beta}_{4}-5r_{h}^{2}\tilde{\beta}_{3}^{2}\right)\right)\Big{]}A_{0}^{\prime
8}$
$\displaystyle-4f^{2}h^{3}r^{2}r_{h}^{6}\Big{[}6h^{2}\bar{\beta}_{4}^{2}\left(13r_{h}^{2}\tilde{\beta}_{3}^{2}-72M_{\text{pl}}^{2}\bar{\beta}_{4}\right)M_{\text{pl}}^{4}+18h^{3}\bar{\beta}_{4}^{2}\left(12M_{\text{pl}}^{2}\bar{\beta}_{4}-5r_{h}^{2}\tilde{\beta}_{3}^{2}\right)M_{\text{pl}}^{4}$
$\displaystyle+\bar{\beta}_{4}\left(3r^{4}\tilde{\beta}_{3}^{2}-8M_{\text{pl}}^{2}\bar{\beta}_{4}\left(8\bar{\beta}_{4}M_{\text{pl}}^{2}+3r_{h}^{2}\tilde{\beta}_{3}^{2}\right)\right)M_{\text{pl}}^{2}$
$\displaystyle+3h\left(\left(r_{h}^{2}\tilde{\beta}_{3}^{4}-M_{\text{pl}}^{2}\bar{\beta}_{4}\tilde{\beta}_{3}^{2}\right)r^{4}+8M_{\text{pl}}^{4}\bar{\beta}_{4}^{2}\left(12\bar{\beta}_{4}M_{\text{pl}}^{2}+r_{h}^{2}\tilde{\beta}_{3}^{2}\right)\right)\Big{]}A_{0}^{\prime
6}$
$\displaystyle+f^{3}h^{2}r^{4}M_{\text{pl}}^{2}r_{h}^{4}\Big{[}-\tilde{\beta}_{3}^{2}r^{4}-288hM_{\text{pl}}^{4}\bar{\beta}_{4}^{2}+24M_{\text{pl}}^{2}\bar{\beta}_{4}\left(4\bar{\beta}_{4}M_{\text{pl}}^{2}+r_{h}^{2}\tilde{\beta}_{3}^{2}\right)$
$\displaystyle+24h^{2}M_{\text{pl}}^{2}\bar{\beta}_{4}\left(9M_{\text{pl}}^{2}\bar{\beta}_{4}-2r_{h}^{2}\tilde{\beta}_{3}^{2}\right)\Big{]}A_{0}^{\prime
4}$
$\displaystyle+2f^{4}hr^{6}M_{\text{pl}}^{4}r_{h}^{2}\left(8\bar{\beta}_{4}M_{\text{pl}}^{2}+r_{h}^{2}\tilde{\beta}_{3}^{2}+h\left(r_{h}^{2}\tilde{\beta}_{3}^{2}-12M_{\text{pl}}^{2}\bar{\beta}_{4}\right)\right)A_{0}^{\prime
2}+f^{5}r^{8}M_{\text{pl}}^{6}\Big{\\}}^{2}$
$\displaystyle\times\Big{\\{}4fhr^{2}A_{0}^{\prime
2}r_{h}^{2}\left[h\left(2\tilde{\beta}_{3}^{2}r_{h}^{2}-3\bar{\beta}_{4}M_{\text{pl}}^{2}\right)+2\bar{\beta}_{4}M_{\text{pl}}^{2}\right]$
$\displaystyle-4h^{2}(3h-2)\bar{\beta}_{4}A_{0}^{\prime
4}r_{h}^{4}\left(h\left(2\tilde{\beta}_{3}^{2}r_{h}^{2}-3\bar{\beta}_{4}M_{\text{pl}}^{2}\right)+2\bar{\beta}_{4}M_{\text{pl}}^{2}\right)+f^{2}r^{4}M_{\text{pl}}^{2}\Big{\\}}^{-3}\,.$

The terms in the first line of Eq. (128) is non-negative as long as
${{\tilde{K}}_{33}}>0$. The positiveness of $c_{\Omega 3{\rm,even}}^{2}$ is
solely determined by the part to the $-3$ power, i.e., the last two lines in
Eq. (128). By treating $f$, $h$, ${\tilde{\beta}}_{3}$, ${\tilde{\beta}}_{4}$
and $A_{0}^{\prime}$ as free parameters, and considering the known constraints
on them, one can determine its minimum by using numerical methods. Taking
advantage of its dimensionless form, this calculation could be fulfilled by
using the built-in functions of Mathematica. This manipulation shows that the
condition $c_{\Omega 3{\rm,even}}^{2}>0$ is always guaranteed.

We proceed to the first and second angular Laplacian stability conditions in
which the story becomes astonishingly sophisticated, in comparing to that of
Model 1. To conquer this problem, we shall first run the analysis around the
$r=r_{h}$ where the algebraic equations in Eq. (126) can be solved for
$A_{0}^{\prime}$ as $A_{0}^{\prime}(r=r_{h})=\sqrt{2\mu}=\kappa$ by choosing
the unit system so that $M_{\text{pl}}=r_{h}$. On using this solution at the
horizon, we start to expand $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ step by
step around the horizon. After a sequence of time-consuming but
straightforward manipulations, their expressions at the horizon converted into
the dimensionless form were obtained. The one for $\mathcal{M}_{1}$ is given
below

$\displaystyle\left.\mathcal{M}_{1}\right|_{r\to r_{h}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{\\{}10\kappa^{12}\tilde{\beta}_{3}^{4}\tilde{\beta}_{4}+4\kappa^{10}\tilde{\beta}_{3}^{2}\tilde{\beta}_{4}\left(2\tilde{\beta}_{3}^{2}+7\tilde{\beta}_{4}\right)+\kappa^{8}\left[\left(12-56\tilde{\beta}_{4}\right)\tilde{\beta}_{3}^{4}+\tilde{\beta}_{4}\left(7-72\tilde{\beta}_{4}\right)\tilde{\beta}_{3}^{2}-16\tilde{\beta}_{4}^{3}\right]$
(129)
$\displaystyle-2\kappa^{6}\left[12\tilde{\beta}_{3}^{4}+3\tilde{\beta}_{4}\left(16\tilde{\beta}_{4}+3\right)\tilde{\beta}_{3}^{2}+4\tilde{\beta}_{4}^{2}\left(24\tilde{\beta}_{4}+1\right)\right]-\kappa^{4}\tilde{\beta}_{4}\left(56\tilde{\beta}_{3}^{2}+160\tilde{\beta}_{4}+1\right)$
$\displaystyle-4\kappa^{2}\left(2\tilde{\beta}_{3}^{2}+11\tilde{\beta}_{4}\right)-4\Big{\\}}\left(4\kappa^{2}\tilde{\beta}_{4}+1\right)^{-2}\left[\kappa^{4}\tilde{\beta}_{3}^{2}-2\kappa^{2}\left(\tilde{\beta}_{3}^{2}+2\tilde{\beta}_{4}\right)-1\right]^{-1}\,,$

while the one for $\mathcal{M}_{2}$ is shown in Appendix B due to its
intricate structure. With the known constraints on ${\tilde{\beta}}_{3}$,
${\tilde{\beta}}_{4}$ and $\mu$ (so does $\kappa$), we can determine the
minimum of $\left.\mathcal{M}_{1}\right|_{r\to r_{h}}$ and
$\left.\mathcal{M}_{2}\right|_{r\to r_{h}}$ on the use of the built-in
functions of Mathematica, which shows that both of them are positive so that
the first as well as second angular Laplacian stability conditions are
satisfied for the given parameter regions.

To further confirm these stability conditions, $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$ will also be plotted out on $r\in[r_{h},+\infty)$ with
chosen coupling constants. This could be fulfilled by choosing ${\mu}=0.5$,
${\tilde{\beta}}_{3}=0.441$ and ${\tilde{\beta}}_{4}=0.1$ as well as being
given the corresponding numerical background solutions shown in the panel (a)
of Fig. 4 in Appendix C. The fundamental steps and techniques of obtaining
these numerical solutions are the same as what we used to draw the panel (a)
of Fig. 3. The only difference is that, for Model 3 there is no analytic
expression for $A_{0}^{\prime}$. Thus, $f$, $h$, and $A_{0}^{\prime}$ need to
be numerically solved simultaneously. For this reason, all of these three
variables are plotted in the panel (a) of Fig. 4. Again, to show a larger
scope, all the quantities are plotted as functions of $\xi$. We note that
$A_{0}^{\prime}$ is plotted in its dimensionless form $r_{h}^{-1}dA_{0}/d\xi$,
which, as shown in the plot, equals $-\xi^{-2}A_{0}^{\prime}$. It’s worth
mentioning here that, $A_{0}^{\prime}$ behaves as ${\cal O}(r^{-2})$ at the
spatial infinity ($\xi\to 0$). That’s why the quantity
$-\xi^{-2}A_{0}^{\prime}$ won’t blow up as $\xi\to 0$. The results of
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are exhibited in the panel (b) of Fig.
4. It is clear that, both these two quantities stay positive at any distances
outside the horizon. As a result, we conclude that this figure further
supports the first and second angular Laplacian stability conditions in Model
3 inside the given parameter space.

|  |  |  Model 1: |  |  Model 2: |  |  Model 3: |
---|---|---|---|---|---|---|---|---
Category |  $\\#$ |  Described by | Type |  See | Type |  See | Type |  See
| first | (70) | Type I | N/A | Type I | N/A | Type I | N/A
No-ghost | second | (72) | Type I | N/A | Type I | N/A | Type I | N/A
| third | (75) | Type II | (114) | Type II | (119) | Type II | (125)
Radial | first | (79) | Type I | N/A | Type I | N/A | Type I | N/A
Laplacian | second | (80) | Type II | (116) | Type II | (121) | Type II | (127)
| third | (78) | Type I | N/A | Type I | N/A | Type I | N/A
Angular | third | (91) | Type II | (117) | Type II | (LABEL:aLap1M2) | Type III | (128)
Laplacian | first | (93) | Type III | (132) | Type IV |  (123) & Fig.3 | Type IV |  (129) & Fig.4
| second | (94) | Type III |  (118) & Fig.1 | Type IV | (134) & Fig.3 | Type IV | (135) & Fig.4
Table 1:  Summarize the stability analysis and corresponding conclusions for
different categories of stability conditions as well as different models. Type
I: This stability condition will automatically get satisfied, provided the
other stability conditions (including those from the odd-parity sector); Type
II: This stability condition can be determined immediately by observing the
corresponding mathematical structure; Type III: The analytic investigation
will be quite complicated so that semi-analytic or numerical techniques have
to be applied in confirming this stability condition. Notice that, the
constraints on coupling parameters, i.e.,
${\tilde{\beta}_{3}}\in(-0.441,0.441)$ and
${\tilde{\beta}_{4}}\in(-0.125,0.251)$ will be considered; Type IV: A complete
numerical analysis is running out of the computational resources. Thus, the
analysis will be carried out by doing expansions around the event horizon
($r=r_{h}$) and also by plotting out those key quantities for certain chosen
coupling parameters (under the constraints mentioned above).

To summary, here we have dealt with 27 cases in total (combining Secs. IV and
V), i.e., 9 stability conditions (3 no-ghost conditions + 3 radial Laplacian
stability conditions + 3 angular Laplacian stability conditions) times 3
distinct models. Basically, no instability was recognized up to now, given the
known$\char 92\relax$updated constraints
${\tilde{\beta}_{3}}\in(-0.441,0.441)$ and
${\tilde{\beta}_{4}}\in(-0.125,0.251)$. For the convenience of readers, we
summarize those main equations, results, analysis, etc. for these 27 cases in
Table 1. Notice that, the main results for each case are categorized in four
different types, denoted by Type I, II, III and IV, respectively, mainly
according to the explicitness of the stability in each case. Type I means the
corresponding stability condition automatically get satisfied, at least by
considering the known constraints to the theory; Type II means the
corresponding stability condition deserves a little analysis. Nonetheless, the
stability is quite explicit so that we can tell that just by looking at the
mathematical structure of the corresponding expression; Type III is a little
sophisticated, so that pure analytic investigation will not work. Even in such
case, we could still confirm these stability conditions by semi-analytic or
numerical methods and with the assistance of computers; Type IV is a little
vague from certain point of view. The problems are really complicated so we
have to expediently check the corresponding stability conditions mainly by
plotting out the results for chosen coupling parameters. A more complete
analysis may be executed in the future for this type of problems.

##  VI Conclusions

In this paper, the stability of static and spherically symmetric BHs with
nontrivial scalar hair in the $U(1)$ GI SVT theory (1) is studied. For the
background ansatz (8), (9), and (10), we considered the even-parity
perturbations given by Eqs. (39), (40), and (41). As a successor of the
previous work Heisenberg  _et al._ (2018), we investigated three types of
stability conditions, i.e., the no-ghost, radial Laplacian and angular
Laplacian ones. In addition to the general analysis, some of the stability
conditions are discussed by considering three distinct concrete models and by
applying suitable analytic, semi-analytic as well as numerical techniques.

To make it more complete, the background field equations (12)-(15) and the
stability conditions against the odd-parity perturbations are revisited in the
presence of the quantities $Y$ and $\tilde{F}$ which are omitted in the
previous work.

On top of that, we first run the general analysis. The original Lagrangian is
expanded to second order for the even-parity perturbations as Eq. (49), which
contains 7 apparent DOFs. After that, with the help of integration by parts,
the integral properties of spherical harmonics, and by introducing suitable
new variables [given in Eqs. (53) and (62)], all the non-dynamical terms
disappear and we are left with a dramatically simplified Lagrangian (63),
adopting 3 DOFs as expectedHeisenberg  _et al._ (2018). On using the resultant
Lagrangian, we derived 9 stability conditions, i.e., 3 no-ghost conditions
(67)-(69) \+ 3 radial Laplacian stability conditions (81) + 3 angular
Laplacian stability conditions (91), (93)-(94).

It is then found that the first [cf., (70)] and second [cf., (72)] no-ghost
conditions as well as the first [cf., (79)] and third [cf., (78)] radial
Laplacian stability conditions get satisfied automatically, provided the other
stability conditions, including those from the odd-parity sector. We then move
to three specific concrete models characterized by Eqs. (109)-(111)] to run
the analysis in Sec. V. By inserting the corresponding functions of each model
into the mathematical expression of the third no-ghost condition, it is found
that this condition hold in all the three models [cf., (114), (119) and
(125)]. Similarly, the second radial Laplacian condition also hold for all
these three models [cf., (116), (121) and (127)].

The situation becomes a little complicated when moving to the angular
Laplacian conditions. First of all, we can easily confirm the third angular
Laplacian condition in an analytic way for Model 1 and Model 2 as in Eqs.
(117) and (LABEL:aLap1M2), respectively. However, for the first and second
angular Laplacian conditions of Model 1, semi-analytic as well as numerical
techniques become necessary. It is found that these stability conditions are
guaranteed by applying the constraints to the coupling parameters obtained
from the odd-parity sector in Ref. Heisenberg  _et al._ (2018) together with
the additional constraints indicated in Fig. 1 (see Appendix C), namely,
${\tilde{\beta}_{3}}\in(-0.441,0.441)$ and
${\tilde{\beta}_{4}}\in(-0.125,0.251)$. With these constraints, the third
angular Laplacian condition for Model 3 can also be confirmed [cf., (128)].
Finally, we notice that the first and second angular Laplacian conditions for
Model 2 and Model 3 are quite difficult to investigate (especially for the
Model 3). To conquer this difficulty, by mimicking Ref. Heisenberg  _et al._
(2018), specific values for the coupling parameters are chosen for these two
models. As a result, we are able to plot out the discriminants [cf., (93) and
(94)] for each case. These two stability conditions are then confirmed in
Figs. 3 and 4 (cf., Appendix C). The main formulas and analysis about these
stability conditions are summarized for the three models in Table 1.

Our current work can be enlarged to several other directions. For instance, we
can dig further into the angular Laplacian conditions in order to put more
stringent constraints on the phase space of the coupling parameters as well as
consider more types of models for analysis. It is also of interest to put
constraints on the model parameters through the observational data of BH
shadow given by the EHT and the other observations. The constraints on the
charges of the supermassive compact objects being black hole candidates are
studied in Refs. Akiyama  _et al._ (2019a, b) assuming that the observed rings
correspond to the photon sphere, as well as in Ref. Tsukamoto and Kase (2024)
assuming that they correspond to the lensing ring. These procedures enable one
to put constraints on the model parameter at the background level. On the
other hand, the second-order Lagrangian analysis is closely related to the
quasi-normal mode problems (see, e.g., Ref. Berti  _et al._ (2018a)). We can
also try to do some relative calculations based on Eq. (63) and put the
further constraint on the model parameters at the perturbation level.

## Acknowledgements

The authors are deeply grateful to Shinji Tsujikawa for useful discussion. CZ
is supported by the National Natural Science Foundation of China under Grant
No. 12205254 and the Science Foundation of China University of Petroleum,
Beijing under Grant No. 2462024BJRC005. RK is supported by the Grant-in-Aid
for Early-Career Scientists of the JSPS No. 20K14471 and the Grant-in-Aid for
Scientific Research (C) of the JSPS No. 23K034210.

##  Appendix A Coefficients in the second-order action of even-parity
perturbations

The coefficients in Eqs. (50) and (51) are given by

$\displaystyle a_{1}\equiv 0\,,\qquad
a_{2}=\frac{r(f^{\prime}h-fh^{\prime})}{fh\phi^{\prime}}a_{4}+\frac{A_{0}^{\prime}}{2}v_{4}+\frac{rA_{0}^{\prime}}{2\phi^{\prime}}v_{6}\,,\qquad
a_{3}=-ra_{4}\,,\qquad a_{4}=\frac{\sqrt{fh}M_{\rm
pl}^{2}}{2}=2f\alpha_{1}\,,$ $\displaystyle
a_{5}=a_{2}^{\prime}-\frac{A_{0}^{\prime\prime}}{2}v_{4}+\frac{A_{0}^{\prime}}{2}(v_{5}-v_{4}^{\prime})\,,\qquad
a_{6}=\frac{r(f^{\prime}h+fh^{\prime})-2fh}{2rfh^{2}\phi^{\prime}}a_{4}-\frac{1}{h\phi^{\prime}}a_{4}^{\prime}-\frac{A_{0}^{\prime}}{4h\phi^{\prime}}v_{6}+\frac{2rf}{\phi^{\prime}}\alpha_{7}\,,$
$\displaystyle
a_{7}=-(ra_{4})^{\prime}-\frac{A_{0}^{\prime}}{4}(2A_{0}^{\prime}v_{1}+\phi^{\prime}v_{4})-\frac{rA_{0}^{\prime}}{2}v_{6}\,,\qquad
a_{8}=-\frac{a_{4}}{2h}\,,\qquad
a_{9}=a_{4}^{\prime}-\frac{rf^{\prime}-2f}{2rf}a_{4}+\frac{A_{0}^{\prime}}{4}v_{6}\,,$
$\displaystyle b_{1}=\frac{a_{4}}{2f}\,,\qquad b_{2}=0\,,\qquad
b_{3}=-\frac{2}{f}a_{2}+\frac{A_{0}^{\prime}}{f}v_{4}\,,\qquad
b_{4}=\frac{2ra_{4}}{f}\,,\qquad b_{5}=-\frac{a_{4}}{f}\,,\qquad c_{1}\equiv
0\,,$ $\displaystyle
c_{2}\equiv\frac{r^{2}\sqrt{fh}\phi^{\prime}(f_{2,X}-h\phi^{\prime
2}f_{2,XX})}{2}-\frac{r^{2}h^{5/2}\phi^{\prime}A_{0}^{\prime 4}(4h\phi^{\prime
2}f_{2,YY}-f_{2,FY})}{f^{3/2}}-\frac{h^{3/2}A_{0}^{\prime
2}}{2\sqrt{f}}\times$ $\displaystyle\hskip 21.33955pt\Big{[}h^{3}\phi^{\prime
5}f_{4,XXX}+2rh^{2}\phi^{\prime 4}f_{3,XX}+h\phi^{\prime
3}\\{6r^{2}f_{2,XY}+(4-13h)f_{4,XX}\\}-12rh\phi^{\prime 2}f_{3,X}$
$\displaystyle\hskip
21.33955pt-\phi^{\prime}\\{r^{2}(6f_{2,Y}+f_{2,XF})+6(2-5h)f_{4,X}-20h\tilde{f}_{4}\\}+6rf_{3}\Big{]}\,,$
$\displaystyle c_{3}=-\frac{1}{2\sqrt{fh}}\frac{\partial{\cal
E}_{11}}{\partial\phi}\,,\qquad c_{4}\equiv-\frac{\sqrt{h}A_{0}^{\prime
2}}{2r\sqrt{f}}\left[h^{2}\phi^{\prime 3}f_{4,XX}+rh\phi^{\prime
2}f_{3,X}-3h\phi^{\prime}(f_{4,X}+2\tilde{f}_{4})-rf_{3}\right]\,,$
$\displaystyle
c_{5}=-\frac{f^{\prime}}{2f}a_{4}-\frac{\phi^{\prime}}{2}d_{2}-\frac{A_{0}^{\prime}}{2}v_{6}+2r\alpha_{6}\,,\qquad
c_{6}=\frac{rf^{\prime}}{2f}a_{4}-\frac{\phi^{\prime}}{4}c_{2}+\frac{r\phi^{\prime}}{4}d_{2}+\frac{A_{0}^{\prime
2}}{4}v_{1}+\frac{\phi^{\prime}A_{0}^{\prime}}{8}v_{4}+\frac{rA_{0}^{\prime}}{2}v_{6}-r^{2}\alpha_{6}\,,$
$\displaystyle d_{1}=\frac{a_{4}}{2f}\,,\qquad
d_{2}\equiv-\frac{h^{3/2}A_{0}^{\prime 2}}{r\sqrt{f}}\left[h^{2}\phi^{\prime
3}f_{4,XX}+rh\phi^{\prime
2}f_{3,X}-2h\phi^{\prime}(3f_{4,X}+2\tilde{f}_{4})-rf_{3}\right]\,,$
$\displaystyle
d_{3}=\frac{2(f^{\prime}h-fh^{\prime})}{rfh\phi^{\prime}}a_{4}+\frac{A_{0}^{\prime}}{r\phi^{\prime}}v_{6}+\frac{A_{0}^{\prime}}{2}\frac{\partial
v_{6}}{\partial\phi}+\frac{4}{h\phi^{\prime}}\alpha_{6}+\frac{4f}{\phi^{\prime}}\alpha_{7}\,,\qquad
d_{4}=-2\alpha_{6}\,,$ $\displaystyle
e_{1}=\frac{1}{fh\phi^{\prime}}\left(a_{2}-2rha_{6}-\frac{A_{0}^{\prime}}{2}v_{4}\right)\,,\qquad
e_{2}=-\frac{1}{\phi^{\prime}}\left(c_{2}+rd_{2}+\frac{A_{0}^{\prime}}{2}v_{4}\right)\,,\qquad
e_{3}=\frac{1}{2}\frac{\partial{\cal E}_{\phi}}{\partial\phi}\,,$
$\displaystyle e_{4}=\frac{1}{r^{2}h\phi^{\prime
2}}\left(1-\frac{3rf^{\prime}}{2f}+\frac{rh^{\prime}}{2h}\right)a_{4}+\frac{h\phi^{\prime\prime}+h^{\prime}\phi^{\prime}}{h\phi^{\prime
2}}c_{4}+\frac{c_{4}^{\prime}}{\phi^{\prime}}+\frac{fA_{0}^{\prime\prime}-2f^{\prime}A_{0}^{\prime}}{4fh\phi^{\prime
2}}v_{6}+\frac{A_{0}^{\prime}}{4h\phi^{\prime 2}}v_{6}^{\prime}$
$\displaystyle\hskip 21.33955pt-\frac{A_{0}^{\prime 2}}{\phi^{\prime
2}}v_{9}-\frac{A_{0}^{\prime 2}}{fh\phi^{\prime
2}}v_{10}-\frac{A_{0}^{\prime}}{\phi^{\prime}}v_{13}-\frac{2rh^{\prime}+2(1-3h)}{h^{2}\phi^{\prime
2}}\alpha_{6}+\frac{2r}{h\phi^{\prime
2}}\alpha_{6}^{\prime}+\frac{2(rf^{\prime}h-f)}{h\phi^{\prime
2}}\alpha_{7}\,,$ $\displaystyle v_{1}\equiv\frac{r^{2}h^{3/2}A_{0}^{\prime
2}}{2f^{3/2}}\left[f_{2,FF}-4h\phi^{\prime 2}(f_{2,FY}-h\phi^{\prime
2}f_{2,YY})\right]$ $\displaystyle\hskip
21.33955pt-\frac{\sqrt{h}}{2\sqrt{f}}\left[2h\phi^{\prime
2}\\{r^{2}f_{2,Y}-h(f_{4,X}+2\tilde{f}_{4})\\}-4rh\phi^{\prime}f_{3}-8(1-h)f_{4}-r^{2}f_{2,F}\right]\,,$
$\displaystyle v_{2}=A_{0}^{\prime}v_{1}\,,\qquad
v_{3}=-A_{0}^{\prime}v_{1}-\frac{\phi^{\prime}}{2}v_{4}-rv_{6}\,,$
$\displaystyle v_{4}\equiv\frac{2r^{2}h^{5/2}\phi^{\prime}A_{0}^{\prime
3}(2h\phi^{\prime
2}f_{2,YY}-f_{2,FY})}{f^{3/2}}+\frac{h^{3/2}A_{0}^{\prime}}{\sqrt{f}}\Big{[}2h\phi^{\prime
3}(r^{2}f_{2,XY}-hf_{4,XX})-4rh\phi^{\prime 2}f_{3,X}$ $\displaystyle\hskip
21.33955pt-r^{2}\phi^{\prime}(4f_{2,Y}+f_{2,XF})+4\phi^{\prime}\\{(3h-2)f_{4,X}+2h\tilde{f}_{4}\\}+4rf_{3}\Big{]}\,,$
$\displaystyle
v_{5}\equiv\frac{\sqrt{h}A_{0}^{\prime}}{\sqrt{f}}\left[2h\phi^{\prime
2}\\{h(f_{4,\phi X}+2\tilde{f}_{4,\phi})-r^{2}f_{2,\phi
Y}\\}+4rh\phi^{\prime}f_{3,\phi}+8(1-h)f_{4,\phi}+r^{2}f_{2,\phi F}\right]\,,$
$\displaystyle v_{6}=4\alpha_{2}\,,\qquad v_{7}=\frac{A_{0}^{\prime
2}}{4}v_{1}\,,\qquad
v_{8}=\frac{1}{2rf}\left[(rf^{\prime}-2f)v_{6}+4rA_{0}^{\prime}v_{10}\right]\,,$
$\displaystyle v_{9}=\alpha_{4}\,,\qquad v_{10}=\alpha_{5}\,,\qquad
v_{11}=-\frac{v_{6}}{2}\,,\qquad v_{12}=-\frac{v_{6}}{2h}\,,$ $\displaystyle
v_{13}\equiv-\frac{2h^{3/2}\phi^{\prime\prime}(h\phi^{\prime
2}f_{4,XX}+r\phi^{\prime}f_{3,X}-f_{4,X}-2\tilde{f}_{4})}{r\sqrt{f}}+\frac{2\sqrt{h}A_{0}^{\prime\prime}[h\phi^{\prime}(f_{4,X}+2\tilde{f}_{4})+rf_{3}]}{r\sqrt{f}}$
$\displaystyle\hskip
21.33955pt-\frac{\sqrt{h}f^{\prime}A_{0}^{\prime}[h\phi^{\prime}(f_{4,X}+2\tilde{f}_{4})+rf_{3}]}{rf^{3/2}}-\frac{h^{\prime}A_{0}^{\prime}[h^{2}\phi^{\prime
3}f_{4,XX}+rh\phi^{\prime
2}f_{3,X}-3h\phi^{\prime}(f_{4,X}+2\tilde{f}_{4})-rf_{3}]}{r\sqrt{fh}}$
$\displaystyle\hskip 21.33955pt-\frac{2\sqrt{h}A_{0}^{\prime}[h\phi^{\prime
2}(\tilde{f}_{3}-2\tilde{f}_{4,\phi}-f_{4,\phi
X})+r\phi^{\prime}(f_{2,Y}-f_{3,\phi})-f_{3}]}{r\sqrt{f}}\,,$ (130)

where ${\cal E}_{11}$ and ${\cal E}_{\phi}$ are the background equations given
in Eqs. (12) and (14), respectively, and $\alpha_{i}$ $(i=1,2,4,5,6,7)$ are
the coefficients in the second-order action for the odd-parity perturbations
defined in Eq. (24).

On the other hand, it’s worth mentioning here that, the coefficient $v_{6}$
adopts a useful feature that

$\displaystyle
v_{6}^{\prime}=\frac{a_{4}\left[2f^{2}\left(rh^{\prime}+2h\right)+hr^{2}\left(f^{\prime}\right)^{2}-fr\left(rf^{\prime}h^{\prime}+2h\left(rf^{\prime\prime}+f^{\prime}\right)\right)\right]-fhr^{2}\left[A_{0}^{\prime}\left(4v_{10}A_{0}^{\prime}+v_{6}f^{\prime}\right)+fv_{6}A_{0}^{\prime\prime}+8\alpha_{7}f^{2}\right]}{f^{2}hr^{2}A_{0}^{\prime}}\,.$
(131)

Notice that, for the economy of notations, here we are writing $\bar{\phi}$
and $\bar{A}_{0}$ appearing in Eqs. (40) and (41) simply as $\phi$ and
$A_{0}$, respectively. It should stimulate no confusions since all the above
coefficients are for the perturbation terms and they themselves are,
definitely, of the zeroth order.

##  Appendix B Expressions for certain quantities appearing in Sec.V

The dimensionless form of $\mathcal{M}_{1}$ for Model 1 (cf. Sec.V.1) is given
below

$\displaystyle\mathcal{M}_{1}$ $\displaystyle=$
$\displaystyle\bigg{\\{}2f^{37/2}\mathfrak{Z}^{2}\xi^{6}\left(f(78h+4)-h\mathfrak{Z}^{2}\xi^{2}\right)\tilde{\beta}_{3}^{2}h^{3/2}+56623104f^{15/2}\mathfrak{Z}^{24}\xi^{64}\left(h\mathfrak{Z}^{2}\xi^{2}+8f(4h-1)\right)\tilde{\beta}_{3}^{20}h^{39/2}$
(132)
$\displaystyle+679477248f^{13/2}\mathfrak{Z}^{28}\xi^{72}\tilde{\beta}_{3}^{22}h^{45/2}+1179648f^{17/2}\mathfrak{Z}^{20}\xi^{56}\Big{[}16\left(5\sqrt{fh^{3}}-40\sqrt{fh^{5}}+53\sqrt{fh^{7}}\right)f^{3/2}$
$\displaystyle-43h^{7/2}\mathfrak{Z}^{4}\xi^{4}+32h\mathfrak{Z}^{2}\xi^{2}\sqrt{f}\left(11\sqrt{fh^{5}}-2\sqrt{fh^{3}}\right)\Big{]}\tilde{\beta}_{3}^{18}h^{15}$
$\displaystyle+393216f^{19/2}\mathfrak{Z}^{16}\xi^{48}\Big{[}4\mathfrak{Z}^{2}\xi^{2}\left(11\sqrt{fh^{3}}-100\sqrt{fh^{5}}+233\sqrt{fh^{7}}\right)f^{3/2}+16\left(19h^{3}-30h^{2}+12h-1\right)\sqrt{h}f^{3}$
$\displaystyle-25h^{7/2}\mathfrak{Z}^{6}\xi^{6}+h\mathfrak{Z}^{4}\xi^{4}\sqrt{f}\left(29\sqrt{fh^{3}}-41\sqrt{fh^{5}}\right)\Big{]}\tilde{\beta}_{3}^{16}h^{13}$
$\displaystyle+2048f^{21/2}\mathfrak{Z}^{14}\xi^{42}\Big{[}16\mathfrak{Z}^{2}\xi^{2}\left(-31\sqrt{fh^{3}}+20\sqrt{fh^{5}}+1865\sqrt{fh^{7}}\right)f^{3/2}$
$\displaystyle+256\left(155h^{3}-123h^{2}+24h-2\right)\sqrt{h}f^{3}-57h^{7/2}\mathfrak{Z}^{6}\xi^{6}+16h\mathfrak{Z}^{4}\xi^{4}\sqrt{f}\left(79\sqrt{fh^{3}}-376\sqrt{fh^{5}}\right)\Big{]}\tilde{\beta}_{3}^{14}h^{11}$
$\displaystyle+32f^{25/2}\mathfrak{Z}^{10}\xi^{30}\Big{[}-16f^{2}\left(577h^{2}-64h+17\right)\mathfrak{Z}^{2}\xi^{2}h^{3/2}-128f(32h+5)\mathfrak{Z}^{4}\xi^{4}h^{5/2}+351\mathfrak{Z}^{6}\xi^{6}h^{7/2}$
$\displaystyle+512f^{5/2}\left(\sqrt{fh}+\sqrt{fh^{3}}-78\sqrt{fh^{5}}+295\sqrt{fh^{7}}\right)\Big{]}\tilde{\beta}_{3}^{10}h^{7}$
$\displaystyle+16f^{6}\mathfrak{Z}^{8}\xi^{24}\Big{[}-8f^{19/2}\left(797h^{2}+12h-5\right)\mathfrak{Z}^{2}\xi^{2}h^{3/2}-2f^{17/2}(279h+41)\mathfrak{Z}^{4}\xi^{4}h^{5/2}+27f^{15/2}\mathfrak{Z}^{6}\xi^{6}h^{7/2}$
$\displaystyle+32\left(f^{21/2}\left(1197h^{2}+4h+30\right)h^{3/2}+\sqrt{f^{21}h}\right)\Big{]}\tilde{\beta}_{3}^{8}h^{5}$
$\displaystyle+2f^{5}\mathfrak{Z}^{6}\xi^{18}\Big{[}3f^{19/2}\mathfrak{Z}^{6}\xi^{6}h^{5/2}-4\left(46f^{21/2}h^{3/2}+3\sqrt{f^{21}h}\right)\mathfrak{Z}^{4}\xi^{4}h$
$\displaystyle+64f^{2}\left(f^{21/2}(489h+83)h^{3/2}+6\sqrt{f^{21}h}\right)+4f\mathfrak{Z}^{2}\xi^{2}\left(f^{21/2}(20-981h)h^{3/2}+3\sqrt{f^{21}h}\right)\Big{]}\tilde{\beta}_{3}^{6}h^{4}$
$\displaystyle+512\mathfrak{Z}^{12}\xi^{36}\Big{[}64f^{29/2}\left(802h^{3}-369h^{2}-3h+2\right)h^{19/2}+219f^{23/2}\mathfrak{Z}^{6}\xi^{6}h^{25/2}$
$\displaystyle-32f^{13}\left(19\sqrt{fh^{3}}-90\sqrt{fh^{5}}-298\sqrt{fh^{7}}\right)\mathfrak{Z}^{2}\xi^{2}h^{9}+8\left(4f^{25/2}h^{23/2}-381(fh)^{25/2}\right)\mathfrak{Z}^{4}\xi^{4}\Big{]}\tilde{\beta}_{3}^{12}$
$\displaystyle+\left(4f^{35/2}(2-57h)\mathfrak{Z}^{6}\xi^{14}h^{7/2}-6f^{33/2}\mathfrak{Z}^{8}\xi^{16}h^{9/2}+8f^{8}\left(f^{21/2}(552h+71)h^{3/2}+\sqrt{f^{21}h}\right)\mathfrak{Z}^{4}\xi^{12}h^{2}\right)\tilde{\beta}_{3}^{4}$
$\displaystyle+2\sqrt{f^{41}h}\bigg{\\}}$
$\displaystyle\times\bigg{\\{}339738624f^{17/2}h^{41/2}\xi^{64}\mathfrak{Z}^{24}\tilde{\beta}_{3}^{20}+2f^{37/2}h^{3/2}\xi^{6}\mathfrak{Z}^{2}\tilde{\beta}_{3}^{2}\left(f(38h+2)-h\xi^{2}\mathfrak{Z}^{2}\right)$
$\displaystyle+28311552f^{19/2}h^{35/2}\xi^{56}\mathfrak{Z}^{20}\tilde{\beta}_{3}^{18}\left(4f(h-1)+3h\xi^{2}\mathfrak{Z}^{2}\right)$
$\displaystyle+1024\xi^{36}\mathfrak{Z}^{12}\tilde{\beta}_{3}^{12}\Big{[}16f^{29/2}\left(319h^{2}-158h-17\right)h^{21/2}+16f^{13}h^{11}\xi^{2}\mathfrak{Z}^{2}\left(139\sqrt{fh^{3}}+35\sqrt{fh}\right)$
$\displaystyle-149\xi^{4}\mathfrak{Z}^{4}(fh)^{25/2}\Big{]}+589824f^{21/2}h^{29/2}\xi^{48}\mathfrak{Z}^{16}\tilde{\beta}_{3}^{16}\left(16f^{2}(h-1)^{2}+32fh(4h-1)\xi^{2}\mathfrak{Z}^{2}-5h^{2}\xi^{4}\mathfrak{Z}^{4}\right)$
$\displaystyle+196608f^{23/2}h^{25/2}\xi^{42}\mathfrak{Z}^{14}\tilde{\beta}_{3}^{14}\left(f^{2}\left(52h^{2}-56h+4\right)+fh(107h+13)\xi^{2}\mathfrak{Z}^{2}-12h^{2}\xi^{4}\mathfrak{Z}^{4}\right)$
$\displaystyle+f^{7}h^{2}\xi^{12}\mathfrak{Z}^{4}\tilde{\beta}_{3}^{4}\left(f^{19/2}h^{5/2}\xi^{4}\mathfrak{Z}^{4}-4h\xi^{2}\mathfrak{Z}^{2}\left(37f^{21/2}h^{3/2}+\sqrt{f^{21}h}\right)+4f\left(f^{21/2}(561h+62)h^{3/2}+\sqrt{f^{21}h}\right)\right)$
$\displaystyle+8f^{6}h^{4}\xi^{18}\mathfrak{Z}^{6}\tilde{\beta}_{3}^{6}\left(9f^{19/2}h^{5/2}\xi^{4}\mathfrak{Z}^{4}+8f^{23/2}(527h+78)h^{3/2}+24\sqrt{f^{23}h}-6h\xi^{2}\mathfrak{Z}^{2}\left(83f^{21/2}h^{3/2}+5\sqrt{f^{21}h}\right)\right)$
$\displaystyle+16f^{5}h^{6}\xi^{24}\mathfrak{Z}^{8}\tilde{\beta}_{3}^{8}\Big{[}109f^{19/2}h^{5/2}\xi^{4}\mathfrak{Z}^{4}+48f^{23/2}(381h+38)h^{3/2}+112\sqrt{f^{23}h}$
$\displaystyle-32h\xi^{2}\mathfrak{Z}^{2}\left(87f^{21/2}h^{3/2}+8\sqrt{f^{21}h}\right)\Big{]}$
$\displaystyle+256f^{4}h^{8}\xi^{30}\mathfrak{Z}^{10}\tilde{\beta}_{3}^{10}\Big{[}45f^{19/2}h^{5/2}\xi^{4}\mathfrak{Z}^{4}-32f^{11}\left(-191\sqrt{fh^{5}}+20\sqrt{fh^{3}}+3\sqrt{fh}\right)$
$\displaystyle+4h\xi^{2}\mathfrak{Z}^{2}\left(3\sqrt{f^{21}h}-91f^{21/2}h^{3/2}\right)\Big{]}+\sqrt{f^{41}h}\bigg{\\}}^{-1}\,,$

where

$\displaystyle\mathfrak{Z}(\xi)$ $\displaystyle\equiv$ $\displaystyle
r_{h}^{-1}\frac{dA_{0}}{d\xi}$ (133) $\displaystyle=$
$\displaystyle\frac{\sqrt[3]{6}\left[\sqrt{6}\sqrt{f^{3}h^{9}\left(54h\kappa^{2}\xi^{6}\tilde{\beta}_{3}^{2}+1\right)}-18f^{3/2}h^{5}\kappa\xi^{3}\tilde{\beta}_{3}\right]^{2/3}-6^{2/3}fh^{3}}{12h^{5/2}\xi^{3}\tilde{\beta}_{3}\left[{\sqrt{6}\sqrt{f^{3}h^{9}\left(54h\kappa^{2}\xi^{6}\tilde{\beta}_{3}^{2}+1\right)}-18f^{3/2}h^{5}\kappa\xi^{3}\tilde{\beta}_{3}}\right]^{1/3}}.$

The dimensionless form of $\mathcal{M}_{2}$ at the $r\to r_{h}$ limit for
Model 2 (cf. Sec.V.2) is given below

$\displaystyle\left.\mathcal{M}_{2}\right|_{r\to r_{h}}$ $\displaystyle=$
$\displaystyle\Big{\\{}\left(1-4\tilde{\beta}_{4}\right)^{2}\left(8\tilde{\beta}_{4}+1\right)^{18}-4\kappa^{2}\tilde{\beta}_{4}\left(4\tilde{\beta}_{4}-1\right)\left(128\tilde{\beta}_{4}^{2}+32\tilde{\beta}_{4}-1\right)\left(8\tilde{\beta}_{4}+1\right)^{15}$
(134)
$\displaystyle-64\kappa^{6}\tilde{\beta}_{4}^{3}\left(8\tilde{\beta}_{4}+1\right)^{14}+12\kappa^{4}\tilde{\beta}_{4}^{2}\left(128\tilde{\beta}_{4}^{2}+16\tilde{\beta}_{4}-5\right)\left(8\tilde{\beta}_{4}+1\right)^{14}$
$\displaystyle+64\tilde{\beta}_{4}^{4}\left(8\tilde{\beta}_{4}\kappa+\kappa\right)^{8}\left(8\tilde{\beta}_{4}+1\right)^{4}$
$\displaystyle-2\kappa^{8}\tilde{\beta}_{3}^{6}\left(\kappa^{2}-16\tilde{\beta}_{4}-2\right)^{2}\left(8\tilde{\beta}_{4}+1\right)^{4}$
$\displaystyle\times\left[155648\tilde{\beta}_{4}^{4}-1024\left(8\kappa^{2}-11\right)\tilde{\beta}_{4}^{3}+16\left(7\kappa^{4}+12\kappa^{2}-180\right)\tilde{\beta}_{4}^{2}+8\left(19\kappa^{2}-4\right)\tilde{\beta}_{4}+25\right]$
$\displaystyle+8\kappa^{10}\tilde{\beta}_{3}^{8}\left(\kappa^{2}-16\tilde{\beta}_{4}-2\right)^{2}$
$\displaystyle\times\left[20\tilde{\beta}_{4}\kappa^{4}+\left(-1088\tilde{\beta}_{4}^{2}-64\tilde{\beta}_{4}+9\right)\kappa^{2}+24\left(8\tilde{\beta}_{4}+1\right)^{2}\left(10\tilde{\beta}_{4}-1\right)\right]\left(8\tilde{\beta}_{4}+1\right)^{2}$
$\displaystyle+2\tilde{\beta}_{3}^{2}\Big{[}-2\kappa^{2}\left(4\tilde{\beta}_{4}-1\right)\left(8\tilde{\beta}_{4}+1\right)^{14}+\kappa^{4}\left(3072\tilde{\beta}_{4}^{4}+3840\tilde{\beta}_{4}^{3}-400\tilde{\beta}_{4}^{2}-40\tilde{\beta}_{4}-1\right)\left(8\tilde{\beta}_{4}+1\right)^{10}$
$\displaystyle-64\kappa^{10}\tilde{\beta}_{4}^{3}\left(8\tilde{\beta}_{4}+1\right)^{8}-2\kappa^{6}\tilde{\beta}_{4}\left(4096\tilde{\beta}_{4}^{4}+15360\tilde{\beta}_{4}^{3}+1216\tilde{\beta}_{4}^{2}-120\tilde{\beta}_{4}-5\right)\left(8\tilde{\beta}_{4}+1\right)^{8}$
$\displaystyle+64\kappa^{12}\tilde{\beta}_{4}^{4}\left(8\tilde{\beta}_{4}+1\right)^{6}+4\tilde{\beta}_{4}^{2}\left(8\tilde{\beta}_{4}\kappa+\kappa\right)^{8}\left(256\tilde{\beta}_{4}^{2}+272\tilde{\beta}_{4}-1\right)\Big{]}\left(8\tilde{\beta}_{4}+1\right)^{2}$
$\displaystyle-32\kappa^{12}\tilde{\beta}_{3}^{10}\left(\kappa^{4}-5\left(8\tilde{\beta}_{4}+1\right)\kappa^{2}+6\left(8\tilde{\beta}_{4}+1\right)^{2}\right)^{2}$
$\displaystyle+\tilde{\beta}_{3}^{4}\Big{[}4\kappa^{4}\left(8\tilde{\beta}_{4}+1\right)^{14}+4\kappa^{6}\left(4608\tilde{\beta}_{4}^{3}-1792\tilde{\beta}_{4}^{2}+96\tilde{\beta}_{4}-5\right)\left(8\tilde{\beta}_{4}+1\right)^{10}$