#  Black hole ringdown from physically sensible initial value problem in
higher-order scalar-tensor theories

Keisuke Nakashi  Department of Social Design Engineering, National Institute
of Technology (KOSEN), Kochi College, 200-1 Monobe Otsu, Nankoku, Kochi,
783-8508, Japan  Department of Physics, Rikkyo University, Toshima, Tokyo
171-8501, Japan  Masashi Kimura  Department of Informatics and Electronics,
Daiichi Institute of Technology, Tokyo 110-0005, Japan  Department of Physics,
Rikkyo University, Toshima, Tokyo 171-8501, Japan  Hayato Motohashi  Division
of Liberal Arts, Kogakuin University, 2665-1 Nakano-machi, Hachioji, Tokyo
192-0015, Japan  Kazufumi Takahashi  Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, 606-8502, Kyoto, Japan

###### Abstract

We study odd-parity perturbations about static and spherically symmetric black
hole solutions with a linearly time-dependent scalar field in higher-order
scalar-tensor theories. In particular, we consider stealth Schwarzschild and
stealth Schwarzschild-de Sitter solutions, where the deviation from the
general relativity case is controlled by a single parameter. We find that
complex frequencies of quasinormal modes (QNMs) are given by a simple scaling
of those in general relativity. We also show that there is a degeneracy
between the parameter characterizing the modification from general relativity
and the black hole mass. We then consider a physically sensible initial value
problem by taking into account the fact that the effective metric for the odd-
parity perturbations is in general different from the background metric. We
confirm that damped oscillations appearing at late times are indeed dominated
by the QNMs. Our analysis includes the case where the perturbations are
superluminal, and we demonstrate in this case that the perturbations can
escape from the region inside the horizon for the background metric.

††preprint: RUP-23-19, YITP-23-118

##  I Introduction

Testing gravity has been a central issue in physics. Apart from cosmological
tests of gravity Koyama:2015vza ; Ferreira:2019xrr ; Arai:2022ilw , there have
been an increasing number of gravitational-wave events from binary black hole
mergers, which offer a possibility to test gravity at strong-field/dynamical
regimes. In general relativity (GR), the late-time gravitational wave signal
emitted from binary black hole mergers, known as the ringdown signal, can be
well described by a superposition of quasinormal modes (QNMs) Buonanno:2006ui
. Each QNM is characterized by a specific complex frequency, whose real and
imaginary parts respectively correspond to the frequency of temporal
oscillation and the exponential damping rate. The no-hair theorem of black
holes in (vacuum) GR implies that the QNM frequencies are determined solely by
the mass and angular momentum of the black hole. However, in modified gravity,
black holes can support some nontrivial hair other than the mass and angular
momentum, which would affect the QNM spectrum. In other words, the information
about the underlying gravitational theory would be encoded in the QNM
spectrum.

In contrast to GR where gravity is described solely by the spacetime metric,
modified gravity theories in general involve additional degrees of freedom.
The simplest class of modified gravity is the class of scalar-tensor theories,
where a single scalar field represents the modification of gravity. Starting
with the seminal theory of Brans-Dicke Brans:1961sx , a number of scalar-
tensor theories have been proposed so far. Horndeski theories Horndeski:1974wa
; Deffayet:2011gz ; Kobayashi:2011nu , which form the most general class of
scalar-tensor theories with second-order Euler-Lagrange equations, provide a
unified description of such traditional theories. It should be noted that the
second-order nature of the Euler-Lagrange equations guarantees the absence of
the Ostrogradsky ghost Woodard:2015zca ; Motohashi:2014opa ; Motohashi:2020psc
; Aoki:2020gfv .

Meanwhile, the Horndeski class is not the most general class of ghost-free
scalar-tensor theories. Indeed, even if the Euler-Lagrange equations contain
higher-order derivatives, the problem of Ostrogradsky ghost can be
circumvented by imposing the degeneracy condition Motohashi:2014opa ;
Langlois:2015cwa ; Motohashi:2016ftl ; Klein:2016aiq ; Motohashi:2017eya ;
Motohashi:2018pxg . Extensions of Horndeski theories in this direction are
called degenerate higher-order scalar-tensor (DHOST) theories Langlois:2015cwa
; Crisostomi:2016czh ; BenAchour:2016fzp . Another systematic way to extend
the Horndeski class is to employ the disformal transformation
Bekenstein:1992pj ; Bruneton:2007si ; Bettoni:2013diz  and its generalization
involving higher derivatives of the scalar field Takahashi:2021ttd ;
Takahashi:2023vva . In fact, the disformal transformation maps the Horndeski
class to (a particular subclass of) the DHOST class, while the generalized
disformal transformation yields a larger class of ghost-free theories, which
is called the generalized disformal Horndeski (GDH) class Takahashi:2022mew
.*1*1*1Matter coupling could introduce an Ostrogradsky mode in generalized
disformal Horndeski theories in general, while there exists a nontrivial
subclass where this problem can be avoided Takahashi:2022mew ; Naruko:2022vuh
; Takahashi:2022ctx ; Ikeda:2023ntu . A yet further extension can be obtained
by relaxing the degeneracy condition in such a way that it is satisfied only
under the unitary gauge. Away from the unitary gauge, apparently there is an
Ostrogradsky mode, but it actually satisfies an elliptic differential equation
on a spacelike hypersurface and hence does not propagate. Such a mode is often
called a shadowy mode DeFelice:2018ewo ; DeFelice:2021hps , which itself is
harmless. By allowing for the existence of the shadowy mode, one obtains
U-DHOST DeFelice:2018ewo ; DeFelice:2021hps ; DeFelice:2022xvq  and
generalized disformal unitary-degenerate (GDU) theories Takahashi:2023jro .

An interesting class of solutions in scalar-tensor theories is the so-called
stealth solution, where the metric is the same as in a GR solution but the
scalar field has a nontrivial profile. The stealth solutions have been found
and studied in the Brans-Dicke theory Nariai1968 ; OHanlon:1972ysn ;
BARROW1990294 ; Romero1993 ; Kolitch:1994kr ; Johri:1994rw ; Giardino:2022sdv
; Giardino:2023qlu , more general scalar-tensor theories Ayon-Beato:2004nzi ;
Ayon-Beato:2005yoq ; Mukohyama:2005rw ; Robinson:2006ib ; Ayon-Beato:2015qfa ;
Alvarez:2016qky ; Smolic:2017bic ; Franzin:2021yvf , and Horndeski and DHOST
theories Babichev:2013cya ; Kobayashi:2014eva ; Babichev:2016kdt ;
Babichev:2017lmw ; Minamitsuji:2018vuw ; BenAchour:2018dap ; Motohashi:2018wdq
; Motohashi:2019sen ; Minamitsuji:2019shy ; Bernardo:2019yxp ;
Charmousis:2019vnf ; Takahashi:2020hso ; Bernardo:2020ehy ; Gorji:2020bfl . In
particular, the general construction of stealth solutions was developed in
Motohashi:2018wdq ; Takahashi:2020hso  in a covariant manner. The perturbation
theory about stealth black hole solutions has been studied extensively
Babichev:2018uiw ; Takahashi:2019oxz ; deRham:2019gha ; Motohashi:2019ymr ;
Khoury:2020aya ; Tomikawa:2021pca ; Takahashi:2021bml ; Mukohyama:2022skk ;
Khoury:2022zor . It then turned out that perturbations of stealth solutions
are strongly coupled in DHOST theories Babichev:2018uiw ; deRham:2019gha ;
Motohashi:2019ymr ; Takahashi:2021bml , and this problem is expected to
persist in GDH theories. A possible way out of this problem is to consider a
small detuning (i.e., scordatura) of the degeneracy condition
Motohashi:2019ymr .*2*2*2The scordatura term affects the stealth black hole
background, leading to a time-dependent correction. However, the time
dependence is typically very weak and can be negligible at astrophysical
scales Mukohyama:2005rw ; DeFelice:2022qaz . This would introduce an
Ostrogradsky mode in general, but its mass can be pushed above the cutoff of
the theory. Moreover, it is even possible to have the scordatura term in
U-DHOST theories that are intrinsically free of Ostrogradsky ghost
DeFelice:2022xvq . Therefore, DHOST (or GDH) theories supplemented with the
scordatura term would provide a consistent description of stealth solutions.

In the present paper, we perform a time-domain analysis of perturbations about
stealth black hole solutions in DHOST theories. In doing so, the main
difficulty comes from the fact that the effective metric (i.e., the one on
which the perturbations propagate) is in general different from the background
metric which determines the motion of (minimally coupled) matter fields. This
implies that a portion of a hypersurface which is spacelike with respect to
the effective metric can be timelike with respect to the background metric.
Therefore, when matter fields are taken into account, one has to carefully
choose the initial hypersurface so that it is spacelike with respect to both
the effective metric and the background metric. This issue has been addressed
in Nakashi:2022wdg  for the case of monopole perturbations about stealth black
hole solutions in DHOST theories. The aim of the present paper is to extend
the analysis of Nakashi:2022wdg  to odd-parity perturbations.

The rest of this paper is organized as follows. In Sec. II, we explain the
DHOST theories and their stealth black hole solutions. In addition, following
Takahashi:2019oxz , we analyze the odd-parity perturbations about the stealth
black hole solutions to see that one has to introduce a new time coordinate
(called $\tilde{t}$) to recast the master equation for the odd-parity
perturbations in the form of a wave equation. In Sec. III, we discuss the
effective metric, the character of a constant-$\tilde{t}$ hypersurface, and
characteristic curves for the odd-parity perturbations about the stealth
Schwarzschild solutions. We also discuss QNM frequencies in the DHOST theories
and obtain the time evolution of the perturbations employing the physically
sensible formulation of an initial value problem developed in Nakashi:2022wdg
. In particular, we confirm that the numerical waveform exhibits damped
oscillations at late times, which can be well fitted by a superposition of the
QNMs for the DHOST theories. In Sec. IV, we perform a similar analysis for the
stealth Schwarzschild-de Sitter solutions. Finally, we draw our conclusions in
Sec. V. In what follows, we use the geometric units in which $c=G=1$.

##  II gravity theory, Background and Odd-parity perturbations

###  II.1 Gravity theory

The action of the quadratic DHOST theories is given by Langlois:2015cwa

$\displaystyle S=\int{\rm
d}^{4}x\sqrt{-g}\left[F_{0}(\phi,X)+F_{1}(\phi,X)\Box\phi+F_{2}(\phi,X)R+\sum_{I=1}^{5}A_{I}(\phi,X)L_{I}^{(2)}\right],$
(1)

where the coupling functions $F_{0},F_{1},F_{2},$ and $A_{I}$ are functions of
the scalar field $\phi$ and its kinetic term $X=\phi_{\mu}\phi^{\mu}$ and

$\displaystyle\begin{split}&L_{1}^{(2)}=\phi_{\mu\nu}\phi^{\mu\nu},\qquad
L_{2}^{(2)}=(\Box\phi)^{2},\qquad
L_{3}^{(2)}=\phi^{\mu}\phi_{\mu\nu}\phi^{\nu}\Box\phi,\\\
&L_{4}^{(2)}=\phi^{\mu}\phi_{\mu\nu}\phi^{\nu\lambda}\phi_{\lambda},\qquad
L_{5}^{(2)}=(\phi^{\mu}\phi_{\mu\nu}\phi^{\nu})^{2},\end{split}$ (2)

with $\phi_{\mu}=\nabla_{\mu}\phi$ and
$\phi_{\mu\nu}=\nabla_{\mu}\nabla_{\nu}\phi$. For a generic choice of the
coupling functions, the theory described by the action (1) suffers from the
problem of the Ostrogradsky ghost associated with higher derivatives in the
equations of motion. The Ostrogradsky ghost can be removed by imposing the
following degeneracy conditions:

$\displaystyle\begin{split}A_{2}&=-A_{1}\neq-\frac{F_{2}}{X},\\\
A_{4}&=\frac{1}{8(F_{2}-XA_{1})^{2}}\left\\{4F_{2}\left[3(A_{1}-2F_{2X})^{2}-2A_{3}F_{2}\right]-A_{3}X^{2}(16A_{1}F_{2X}+A_{3}F_{2})\right.\\\
&\quad\left.+4X(3A_{1}A_{3}F_{2}+16A_{1}^{2}F_{2X}-16A_{1}F_{2X}^{2}-4A_{1}^{3}+2A_{3}F_{2}F_{2X})\right\\},\\\
A_{5}&=\frac{1}{8(F_{2}-XA_{1})^{2}}(2A_{1}-XA_{3}-4F_{2X})\left[A_{1}(2A_{1}+3XA_{3}-4F_{2X})-4A_{3}F_{2}\right],\end{split}$
(3)

where a subscript $X$ denotes the derivative with respect to $X$. The DHOST
theories described by Eq. (1) with the degeneracy conditions (3) is called
class Ia Langlois:2015cwa ; BenAchour:2016cay , which can be mapped to the
Horndeski theory via disformal transformation. It is known that all the other
classes of quadratic DHOST theories are phenomenologically disfavored in the
sense that either the cosmological perturbations are unstable or the modes
correspond to gravitational waves are absent.

In the present paper, we consider a subclass of the class Ia quadratic DHOST
theories, which is described by the following action:

$\displaystyle S=\int{\rm
d}x^{4}\sqrt{-g}\left[F_{0}(X)+F_{2}(X)R+\sum_{I=1}^{5}A_{I}(X)L_{I}^{(2)}\right],$
(4)

where we have set $F_{1}=0$ and assumed that the coupling functions are
functions only of $X$. In other words, we focus on the subclass of the
quadratic DHOST theories whose action is invariant under the shift
($\phi\to\phi+{\rm const}.$) and the reflection ($\phi\to-\phi$) of the scalar
field. As we will see in the next subsection, these theories admit an
interesting class of solutions known as the stealth solutions, i.e., a GR
solution with a linearly time-dependent scalar field.

###  II.2 Background spacetime and scalar field

We consider a static and spherically symmetric background spacetime. The
metric of the background spacetime is given by

$\displaystyle\bar{g}_{\mu\nu}{\rm d}x^{\mu}{\rm d}x^{\nu}=-A(r){\rm
d}t^{2}+\frac{{\rm d}r^{2}}{B(r)}+r^{2}\gamma_{ab}{\rm d}x^{a}{\rm d}x^{b},$
(5)

where $\gamma_{ab}$ is the metric on a two-dimensional unit sphere,
$\gamma_{ab}{\rm d}x^{a}{\rm d}x^{b}={\rm d}\theta^{2}+\sin^{2}\theta{\rm
d}\varphi^{2}$. As for the scalar field, we impose the following ansatz:

$\displaystyle\bar{\phi}(t,r)=qt+\psi(r),$ (6)

where $q$ is a nonvanishing constant. We note that the linear time dependence
of the scalar field is compatible with the static metric because the action
(4) depends on the scalar field only through its derivatives. Having said
that, the linear time dependence can be still allowed in theories without
shift symmetry Minamitsuji:2018vuw ; Motohashi:2018wdq ; Takahashi:2020hso .

In the present paper, in particular, we focus on stealth black hole solutions.
A stealth black hole solution is described by the metric which is the same as
the one in GR, while the scalar field has a nontrivial configuration. The
general construction of stealth solutions was developed in Motohashi:2018wdq ;
Takahashi:2020hso  in a covariant manner. The idea is to substitute the metric
and scalar field ansatz into the equations of motion and derive the conditions
on the coupling functions of DHOST theories under which the equations are
trivially satisfied. Assuming that $X=-q^{2}$, the stealth Schwarzschild-de
Sitter (dS) metric,

$\displaystyle A(r)=B(r)=1-\frac{r_{\rm s}}{r}-\frac{\Lambda r^{2}}{3},$ (7)

with $r_{\rm s}$ and $\Lambda$ being constants, can be a solution if the
following conditions are satisfied Motohashi:2018wdq ; Takahashi:2020hso :

$\displaystyle\left.\left\\{F_{0}+2\Lambda\left(F_{2}-XA_{1}\right)\right\\}\right|_{X=-q^{2}}=0,\,\,\left.\left\\{2F_{0X}+\Lambda\left(8F_{2X}-2A_{1}+4XA_{1X}+3XA_{3}\right)\right\\}\right|_{X=-q^{2}}=0.$
(8)

Note that, among the three degeneracy conditions in (3), we have used only
$A_{2}=-A_{1}$ in deriving the above conditions. Therefore, the stealth
Schwarzschild-dS solution exists even away from the DHOST theories so long as
$A_{2}=-A_{1}$. Note also that the stealth Schwarzschild solution can be
realized by putting $\Lambda=0$. In this case, the above condition reads

$\displaystyle\left.F_{0}\right|_{X=-q^{2}}=0,\qquad\left.F_{0X}\right|_{X=-q^{2}}=0.$
(9)

For the stealth black hole solutions, the scalar field profile can be obtained
from the condition $X=-q^{2}$ as follows:

$\displaystyle\bar{\phi}=q\left(t\pm\int\frac{\sqrt{1-A(r)}}{A(r)}{\rm
d}r\right).$ (10)

Here, we choose the plus branch so that $\phi$ is regular at the future event
horizon. Indeed, for the plus branch, the behavior of the scalar field near
the future event horizon where $A(r)\simeq 0$ can be approximated as

$\displaystyle\bar{\phi}\simeq q\left(t+\int\frac{{\rm d}r}{A(r)}\right)=qv,$
(11)

where $v$ is the ingoing Eddington-Finkelstein coordinate defined by
$v\coloneqq t+\int A(r)^{-1}{\rm d}r$.

###  II.3 Odd-parity perturbations: quadratic Lagrangian and equation of
motion

We study linear odd-parity perturbations around a static and spherically
symmetric spacetime in DHOST theories. Although we will focus on the stealth
black hole solutions in the subsequent sections, for the time being, we
investigate the perturbations around a general static and spherically
symmetric spacetime described by the metric (5), following the discussion in
Takahashi:2019oxz ; Takahashi:2021bml . To study the odd-parity perturbations,
we define the metric perturbation as $\epsilon h_{\mu\nu}\coloneqq
g_{\mu\nu}-\bar{g}_{\mu\nu}$, where $\epsilon$ is a small parameter. Due to
the spherical symmetry of the background spacetime, it is useful to expand the
odd-parity perturbations in terms of the spherical harmonics $Y_{\ell
m}(\theta,\varphi)$ as follows:

$\displaystyle\begin{split}h_{tt}&=h_{tr}=h_{rr}=0,\\\
h_{ta}&=\sum_{\ell,m}h_{0,\ell m}(t,r)E_{a}{}^{b}\bar{\nabla}_{b}Y_{\ell
m}(\theta,\varphi),\\\ h_{ra}&=\sum_{\ell,m}h_{1,\ell
m}(t,r)E_{a}{}^{b}\bar{\nabla}_{b}Y_{\ell m}(\theta,\varphi),\\\
h_{ab}&=\sum_{\ell,m}h_{2,\ell
m}(t,r)E_{(a}{}^{c}\bar{\nabla}_{b)}\bar{\nabla}_{c}Y_{\ell
m}(\theta,\varphi),\end{split}$ (12)

where $E_{ab}$ is the completely antisymmetric tensor defined on a two-
dimensional unit sphere, and $\bar{\nabla}_{a}$ denotes the covariant
derivative with respect to $\gamma_{ab}$. Due to the symmetry of the
background spacetime, it is sufficient to consider only $m=0$. We note that
the odd-parity perturbations do not have $\ell=0$ mode, and $h_{2}$ vanishes
for $\ell=1$. In what follows, we focus on the modes with $\ell\geq 2$ where
the odd-parity perturbations are dynamical. Also, we do not consider the
perturbation of the scalar field, because it belongs to the even-parity
perturbations.

In order to eliminate an unphysical degree of freedom, we consider an
infinitesimal coordinate transformation: $x^{a}\to x^{a}+\epsilon\xi^{a}$. A
general infinitesimal transformation for the odd-parity modes can be written
as

$\displaystyle\xi^{a}=\sum_{\ell,m}\Xi_{\ell
m}(t,r)E^{ab}\bar{\nabla}_{b}Y_{\ell m}(\theta,\varphi).$ (13)

Then, the gauge transformation law for the perturbation variables is given by

$\displaystyle h_{0}\to h_{0}-\dot{\Xi},\qquad h_{1}\to
h_{1}-\Xi^{\prime}+\frac{2}{r}\Xi,\qquad h_{2}\to h_{2}-2\Xi,$ (14)

where a dot and a prime denote the derivatives with respect to $t$ and $r$,
respectively. For $\ell\geq 2$, we set $h_{2}=0$ to fix the gauge freedom,
which is a complete gauge fixing and hence we can legitimately impose it at
the action level Motohashi:2016prk .

The quadratic Lagrangian can be written in terms of a master variable
$\chi_{\ell}$ as follows Takahashi:2019oxz :

$\displaystyle\frac{2\ell+1}{2\pi}{\cal
L}^{(2)}=\frac{\ell(\ell+1)}{2(\ell-1)(\ell+2)}\sqrt{\frac{B}{A}}\left\\{b_{1}\dot{\chi}_{\ell}^{2}-b_{2}\chi_{\ell}^{\prime
2}+b_{3}\dot{\chi}_{\ell}\chi_{\ell}^{\prime}-\left[\ell(\ell+1)b_{4}+V_{\rm
eff}(r)\right]\chi_{\ell}^{2}\right\\},$ (15)

where

$\displaystyle b_{1}=\frac{r^{2}{\cal FH}^{2}}{A{\cal FG}+B{\cal
J}^{2}},\qquad b_{2}=\frac{r^{2}AB{\cal GH}^{2}}{A{\cal FG}+B{\cal
J}^{2}},\qquad b_{3}=\frac{2r^{2}B{\cal H}^{2}{\cal J}}{A{\cal FG}+B{\cal
J}^{2}},\qquad b_{4}={\cal H},$ (16)

and $V_{\rm eff}(r)$ is given by

$\displaystyle V_{\rm eff}(r)=r^{2}{\cal
H}\left[b_{2}\sqrt{\frac{B}{A}}\left(\frac{1}{r^{2}{\cal
H}}\sqrt{\frac{A}{B}}\right)^{\prime}\,\right]^{\prime}-2{\cal H},$ (17)

with ${\cal F}$, ${\cal G}$, ${\cal H}$, and ${\cal J}$ defined by

$\displaystyle\begin{split}&{\cal
F}=2\left(F_{2}+\frac{q^{2}}{A}A_{1}\right),\qquad{\cal
G}=2\left[F_{2}-\left(\frac{q^{2}}{A}+X\right)A_{1}\right],\\\ &{\cal
H}=2\left(F_{2}-XA_{1}\right),\qquad{\cal
J}=-2q\psi^{\prime}A_{1}.\end{split}$ (18)

The relation between the master variable and the original perturbation
variables can be found in Takahashi:2019oxz . The existence of the cross term
$b_{3}\dot{\chi}_{\ell}\chi_{\ell}^{\prime}$ is the crucial difference from
the case with $q=0$, $\psi^{\prime}=0$, and/or $A_{1}=0$. Indeed, we have
$b_{3}\propto{\cal J}\propto q\psi^{\prime}A_{1}$, and hence the cross term
vanishes if $q\psi^{\prime}A_{1}=0$. However, in the present paper, we do not
consider the case where $q\psi^{\prime}A_{1}=0$ because in this case, the
equation of motion and consequently the evolution of the odd-parity
perturbations are completely the same as those in GR.

Let us proceed with the quadratic Lagrangian (15). We can eliminate the cross
term $b_{3}\dot{\chi}_{\ell}\chi_{\ell}^{\prime}$ by introducing a new
coordinate $\tilde{t}$ as follows:

$\displaystyle\tilde{t}=t+\int\frac{b_{3}}{2b_{2}}{\rm d}r.$ (19)

With this new coordinate, the quadratic Lagrangian becomes

$\displaystyle{\cal L}^{(2)}\propto\tilde{{\cal
L}}=\frac{1}{2}\sqrt{\frac{B}{A}}\left\\{\tilde{b}_{1}(\partial_{\tilde{t}}\chi_{\ell})^{2}-b_{2}\chi_{\ell}^{\prime
2}-\left[\ell(\ell+1)b_{4}+V_{\rm eff}(r)\right]\chi_{\ell}^{2}\right\\},$
(20)

where

$\displaystyle\tilde{b}_{1}=b_{1}+\frac{b_{3}^{2}}{4b_{2}}.$ (21)

Next, we obtain the equation of motion for the odd-parity perturbations.
Varying the quadratic Lagrangian (20) with respect to the master variable
$\chi_{\ell}$, we obtain the equation of motion as

$\displaystyle-\partial_{\tilde{t}}^{2}\chi_{\ell}+\frac{b_{2}}{\tilde{b}_{1}}\chi_{\ell}^{\prime\prime}+\frac{Ab_{2}B^{\prime}+B(2Ab_{2}^{\prime}-b_{2}A^{\prime})}{2AB\tilde{b}_{1}}\chi_{\ell}^{\prime}-\frac{\ell(\ell+1)b_{4}+V_{\rm
eff}}{\tilde{b}_{1}}\chi_{\ell}=0.$ (22)

We introduce a new coordinate $\tilde{x}$ and a new variable $\Psi$ to
transform the above equation into the form of a two-dimensional wave equation:

$\displaystyle\tilde{x}$
$\displaystyle=\int\sqrt{\frac{\tilde{b}_{1}}{b_{2}}}{\rm d}r,$ (23)
$\displaystyle\Psi_{\ell}$ $\displaystyle=\frac{\chi_{\ell}}{F(\tilde{x})},$
(24)

where $F(\tilde{x})$ is given by

$\displaystyle F(\tilde{x})=\left(\frac{A}{B\tilde{b}_{1}b_{2}}\right)^{1/4}.$
(25)

Note that $\tilde{x}$ is a generalization of the tortoise coordinate.
Consequently, the equation of motion becomes

$\displaystyle\left[\frac{\partial^{2}}{\partial\tilde{x}^{2}}-\frac{\partial^{2}}{\partial\tilde{t}^{2}}-V_{\ell}(\tilde{x})\right]\Psi_{\ell}=0,$
(26)

where $V_{\ell}(\tilde{x})$ is the effective potential defined by

$\displaystyle V_{\ell}(\tilde{x})=\frac{\ell(\ell+1)b_{4}+V_{\rm
eff}}{\tilde{b}_{1}}+F\frac{{\rm d}^{2}}{{\rm
d}\tilde{x}^{2}}\left(\frac{1}{F}\right).$ (27)

When we fix the background solution, we can compute the effective potential
$V_{\ell}(\tilde{x})$ from the above formula, and hence we can investigate the
time evolution of the odd-parity perturbations based on the master equation
(26).

It should be noted that one can derive a master equation of the same form even
if we do not impose the degeneracy conditions (3), as clarified in
Tomikawa:2021pca . This is as expected because an extra scalar degree of
freedom belongs to the even-parity perturbations and hence does not affect the
odd-parity sector. As mentioned earlier in Sec. II.2, so long as
$A_{2}=-A_{1}$ is satisfied, the class of higher-order scalar-tensor theories
described by the action (4) allows for the stealth Schwarzschild-dS solution
under the condition (8). Moreover, even when $A_{2}\neq-A_{1}$ (which happens
if we take into account the scordatura term Motohashi:2019ymr ), the deviation
of the background solution from the stealth Schwarzschild-dS profile is
typically very weak and can be negligible at astrophysical scales
Mukohyama:2005rw ; DeFelice:2022qaz . Therefore, it is not necessary to impose
the degeneracy conditions (3) for the study of perturbations about the stealth
Schwarzschild-dS profile. Having said that, for concreteness, we focus on the
stealth Schwarzschild(-dS) solution in the DHOST theories in the subsequent
analyses.

##  III Stealth Schwarzschild solutions

###  III.1 Effective metric

In this section, we consider the stealth Schwarzschild profile as the
background solution. From the diagonalized quadratic Lagrangian (20), we can
find the effective metric on which the odd-parity perturbations propagate. In
what follows, we are interested in the propagation of odd-parity perturbations
in the radial direction, and hence we focus on the first two terms in (20) and
define a two-dimensional effective metric $Z_{IJ}$ ($I,J=\\{\tilde{t},r\\}$)
as

$\displaystyle\tilde{{\cal L}}_{\rm
kin}=\sqrt{\frac{B}{A}}\left[\frac{\tilde{b}_{1}}{2}(\partial_{\tilde{t}}\chi_{\ell})^{2}-\frac{b_{2}}{2}\chi_{\ell}^{\prime
2}\right]\eqqcolon-\frac{1}{2}Z^{IJ}\partial_{I}\chi_{\ell}\partial_{J}\chi_{\ell},$
(28)

where $Z^{IJ}$ is the inverse of $Z_{IJ}$. The component of the effective
metric is given by

$\displaystyle Z_{IJ}{\rm d}x^{I}{\rm
d}x^{J}=\sqrt{\frac{A}{B}}\left[-\frac{1}{\tilde{b}_{1}}{\rm
d}\tilde{t}^{2}+\frac{1}{b_{2}}{\rm d}r^{2}\right].$ (29)

Note that the effective metric is in general different from the background
metric, i.e., $Z_{IJ}{\rm d}x^{I}{\rm d}x^{J}\neq\bar{g}_{IJ}{\rm d}x^{I}{\rm
d}x^{J}$. For the stealth Schwarzschild solutions, $Z_{\tilde{t}\tilde{t}}$
becomes

$\displaystyle Z_{\tilde{t}\tilde{t}}=-\frac{F_{2}(r-r_{\rm
s})-q^{2}A_{1}r_{\rm s}}{2r^{3}(F_{2}+q^{2}A_{1})^{2}}.$ (30)

For the spacetime described by the effective metric $Z_{IJ}$, the vector field
$\partial_{\tilde{t}}$ is a Killing vector field. The Killing horizon is
located at the radius where $Z_{\tilde{t}\tilde{t}}$ changes its sign. From
Eq. (30), the radius of the Killing horizon, denoted by $r_{\rm g}$, can be
read off as

$\displaystyle r_{\rm g}=\left(1+\frac{q^{2}A_{1}}{F_{2}}\right)r_{\rm
s}\eqqcolon(1+\zeta)r_{\rm s}.$ (31)

Since the conditions for no ghost/gradient instabilities are given by
Takahashi:2021bml

$\displaystyle F_{2}>0,\qquad F_{2}+q^{2}A_{1}>0,$ (32)

the Killing horizon $r_{\rm g}$ is positive. Note that these conditions imply
$\zeta>-1$. Note also that $r_{\rm g}>r_{\rm s}$ for $\zeta>0$, while $r_{\rm
g}<r_{\rm s}$ for $\zeta<0$. The two radii coincide with each other for
$q^{2}A_{1}=0$, or equivalently $\zeta=0$.

###  III.2 Characters of a constant-$\tilde{t}$ surface

Next, we discuss characters of the new time coordinate $\tilde{t}$. For the
stealth Schwarzschild solutions, $\tilde{t}$ can be analytically obtained from
Eq. (19) as follows:

$\displaystyle\tilde{t}=t+2\sqrt{\frac{r}{r_{\rm s}}}(r_{\rm s}-r_{\rm
g})-\frac{1}{\sqrt{r_{\rm s}}}\left(r_{\rm
g}^{3/2}\log\left|\frac{\sqrt{r}-\sqrt{r_{\rm g}}}{\sqrt{r}+\sqrt{r_{\rm
g}}}\right|-r_{\rm s}^{3/2}\log\left|\frac{\sqrt{r}-\sqrt{r_{\rm
s}}}{\sqrt{r}+\sqrt{r_{\rm s}}}\right|\right)+\tilde{t}_{\rm c},$ (33)

where $\tilde{t}_{\rm c}$ is an integration constant. Let us investigate
whether a constant-$\tilde{t}$ surface is spacelike with respect to the
background metric or not. To this end, we consider a vector field
$\partial_{\mu}\tilde{t}$ which is normal to a constant-$\tilde{t}$ surface.
The norm of $\partial_{\mu}\tilde{t}$ associated with the background metric is
given by

$\displaystyle\bar{g}^{\mu\nu}\partial_{\mu}\tilde{t}\partial_{\nu}\tilde{t}=\frac{r(r_{\rm
g}^{2}-rr_{\rm s})}{r_{\rm s}(r-r_{\rm g})^{2}}.$ (34)

Therefore, the constant-$\tilde{t}$ surface is spacelike for $r>r_{\rm
g}^{2}/r_{\rm s}$, while it is timelike for $r<r_{\rm g}^{2}/r_{\rm s}$. Now,
we discuss the relation between the location of the Killing horizon for the
odd-parity perturbations $r_{\rm g}$ and the characteristic radius $r_{\rm
g}^{2}/r_{s}$. The Killing horizon $r_{\rm g}$ is greater than the
characteristic radius $r_{\rm g}^{2}/r_{\rm s}$ if $r_{\rm s}>r_{\rm g}$, or
equivalently $\zeta<0$. Consequently, if we focus on the spacetime in the
range $r>r_{\rm g}$, the constant-$\tilde{t}$ surface is always spacelike. On
the other hand, the Killing horizon $r_{\rm g}$ is smaller than the
characteristic radius $r_{\rm g}^{2}/r_{\rm s}$ if $r_{\rm s}<r_{\rm g}$, or
equivalently $\zeta>0$. Therefore, the constant-$\tilde{t}$ surface becomes
spacelike in the range $r>r_{\rm g}^{2}/r_{\rm s}$, while it becomes timelike
in the range $r_{\rm g}<r<r_{\rm g}^{2}/r_{\rm s}$. Figure 1 shows the typical
behavior of the constant-$\tilde{t}$ surface embedded in the Penrose diagram
of the Schwarzschild spacetime. The black solid curves are the
constant-$\tilde{t}$ surfaces. In the yellow shaded region, the
constant-$\tilde{t}$ surfaces are spacelike.

Figure 1: Typical behavior of constant-$\tilde{t}$ surface for (A) $\zeta>0$
and (B) $\zeta<0$ embedded in the Penrose diagram of the Schwarzschild
spacetime. The black curves represent constant-$\tilde{t}$ surfaces. The
constant-$\tilde{t}$ surface is spacelike in the yellow shaded region.

###  III.3 Characteristic curves

In the high-frequency regime, the odd-parity perturbations propagate along the
characteristic curves on which either $\tilde{v}=\tilde{t}+\tilde{x}={\rm
const}.$ or $\tilde{u}=\tilde{t}-\tilde{x}={\rm const}.$ is satisfied. To
understand properties of the characteristic curves, we perform a similar
analysis as the one in the previous subsection. That is, we study the vector
fields $\partial_{\mu}\tilde{u}$ and $\partial_{\mu}\tilde{v}$ which are
normal to the characteristic curves. The norms of these vector fields with
respect to the background metric are given by

$\displaystyle\bar{g}^{\mu\nu}\partial_{\mu}\tilde{u}\partial_{\nu}\tilde{u}=\frac{r(r_{\rm
g}-r_{\rm s})}{r_{\rm s}(\sqrt{r}-\sqrt{r_{\rm
g}})^{2}},\qquad\bar{g}^{\mu\nu}\partial_{\mu}\tilde{v}\partial_{\nu}\tilde{v}=\frac{r(r_{\rm
g}-r_{\rm s})}{r_{\rm s}(\sqrt{r}+\sqrt{r_{\rm g}})^{2}},$ (35)

respectively. Therefore, for $r_{\rm g}>r_{\rm s}$ or equivalently $\zeta>0$,
the characteristic curves are timelike, while for $r_{\rm g}<r_{\rm s}$ or
equivalently $\zeta<0$, the characteristic curves are spacelike, i.e., the
odd-parity perturbations become superluminal. For $\zeta<0$ case, due to the
superluminal propagation, perturbations can propagate from the region in
$r_{\rm g}<r<r_{\rm s}$ to that in $r>r_{\rm s}$ (see Appendix A).

Figure 2: The characteristic curves for (A) $\zeta>0$ and (B) $\zeta<0$
embedded in the Penrose diagram of the Schwarzschild spacetime. The red curves
and the blues curves represent constant-$\tilde{v}$ curves and
constant-$\tilde{u}$ curves, respectively. For (A) $\zeta>0$, the
characteristic curves are always timelike, while for (B) $\zeta<0$, the
characteristic curves are spacelike, i.e., the odd-parity perturbations are
superluminal.

Figure 2 shows the characteristic curves embedded in the Penrose diagram of
the Schwarzschild spacetime.

###  III.4 Equation of motion and QNM frequencies

Let us study the master equation (26) for the case of stealth Schwarzschild
solutions. The generalized tortoise coordinate $\tilde{x}$ and the new master
variable $\Psi_{\ell}$ defined in Eqs. (23) and (24) take the form of

$\displaystyle\tilde{x}$ $\displaystyle=\sqrt{1+\zeta}\left[r+r_{\rm
g}\log\left|\frac{r}{r_{\rm g}}-1\right|\right],$ (36)
$\displaystyle\Psi_{\ell}$ $\displaystyle=r\sqrt{2F_{2}}\left(\frac{r_{\rm
g}}{r_{\rm s}}\right)^{3/4}\chi_{\ell}.$ (37)

We note that $\tilde{x}\to-\infty$ as $r\to r_{\rm g}$ and
$\tilde{x}\to\infty$ as $r\to\infty$. Here, we have chosen the integration
constant for $\tilde{x}$ so that $\tilde{x}=0$ at $r=0$. The master equation
(26) is now written as

$\displaystyle\left[\frac{\partial^{2}}{\partial\tilde{x}^{2}}-\frac{\partial^{2}}{\partial\tilde{t}^{2}}-V_{\ell}(\tilde{x})\right]\Psi_{\ell}=0,$
(38)

where

$\displaystyle V_{\ell}(\tilde{x})=\frac{1}{1+\zeta}\left(1-\frac{r_{\rm
g}}{r}\right)\left[\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{\rm
g}}{r^{3}}\right].$ (39)

Note that if $\zeta=0$, the above equation reduces to the standard Regge-
Wheeler equation in GR.

It should be noted that the master equation (38) for the odd-parity
perturbations about stealth solutions in the DHOST theory is the same as the
one in GR except that the effective potential is multiplied by the factor of
$(1+\zeta)^{-1}$ [see Eq. (39)]. Indeed, if we introduce rescaled coordinates
$\tilde{T}$ and $\tilde{X}$ as

$\displaystyle\tilde{T}$ $\displaystyle=\frac{\tilde{t}}{\sqrt{1+\zeta}},$
(40) $\displaystyle\tilde{X}$
$\displaystyle=\frac{\tilde{x}}{\sqrt{1+\zeta}}=r+r_{\rm
g}\log\left|\frac{r}{r_{\rm g}}-1\right|,$ (41)

then the master equation (38) can be rewritten as

$\displaystyle\left[\frac{\partial^{2}}{\partial\tilde{X}^{2}}-\frac{\partial^{2}}{\partial\tilde{T}^{2}}-\tilde{V}_{\ell}(\tilde{X})\right]\Psi_{\ell}=0,$
(42)

with

$\displaystyle\tilde{V}_{\ell}(\tilde{X})=\left(1-\frac{r_{\rm
g}}{r}\right)\left[\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{\rm
g}}{r^{3}}\right].$ (43)

Equation (42) is nothing but the standard Regge-Wheeler equation in GR if we
identify $r_{\rm g}$ as the Schwarzschild radius. This implies that we can map
a solution for the wave equation in GR to a solution in the DHOST theory: The
latter is obtained by just rescaling the coordinates in the former. This fact
can be used to discuss the QNM frequencies and the power-law tail in the DHOST
theory.

Let us first discuss the QNM frequencies. Substituting the ansatz
$\Psi_{\ell}=\psi_{\ell}(\tilde{X})e^{-i\tilde{W}\tilde{T}}$ into Eq. (42), we
have

$\displaystyle\left[-\frac{{\rm d}^{2}}{{\rm
d}\tilde{X}^{2}}+\tilde{V}_{\ell}(\tilde{X})\right]\psi_{\ell}(\tilde{X})=\tilde{W}^{2}\psi_{\ell}(\tilde{X}).$
(44)

The QNMs are defined as the modes that are purely ingoing ($\psi_{\ell}\sim
e^{-i\tilde{W}\tilde{X}}$) as $r\to r_{\rm g}$, and purely outgoing
($\psi_{\ell}\sim e^{i\tilde{W}\tilde{X}}$) as $r\to\infty$. Let
$\omega^{\text{Sch}}_{\ell,n}(r_{\rm s})$ be the QNM frequencies for the
Schwarzschild spacetime in GR obtained by solving the standard Regge-Wheeler
equation, where $n$ is the overtone number. For instance,
$\omega_{2,0}^{\text{Sch}}=(0.74734-0.17792\,i)/r_{\rm s}$ for the $\ell=2$
fundamental mode. Also, let $\tilde{W}_{\ell,n}(r_{\rm g})$ be the QNM
frequencies obtained by solving Eq. (44). The relation between
$\omega^{\text{Sch}}_{\ell,n}$ and $\tilde{W}_{\ell,n}$ is given by $r_{\rm
g}\,\tilde{W}_{\ell,n}=r_{\rm s}\,\omega^{\text{Sch}}_{\ell,n}$. From Eq.
(40), we can rewrite the ansatz for $\Psi_{\ell}$ as
$\Psi_{\ell}=\psi_{\ell}(\tilde{x})e^{-\tilde{W}\tilde{t}/\sqrt{1+\zeta}}\eqqcolon\psi_{\ell}(\tilde{x})e^{-i\omega^{\rm
DHOST}\tilde{t}}$. Then, the QNM frequencies $\omega^{\rm DHOST}_{\ell,n}$ can
be expressed in terms of $\omega^{\text{Sch}}_{\ell,n}$ as

$\displaystyle\omega^{\rm DHOST}_{\ell,n}=\frac{r_{\rm
s}\,\omega^{\text{Sch}}_{\ell,n}}{r_{\rm s}(1+\zeta)^{3/2}},$ (45)

where we have used $r_{\rm g}=(1+\zeta)r_{\rm s}$. Note that the numerator is
the QNM frequencies of the Schwarzschild spacetime in GR in unit of $r_{\rm
s}^{-1}$, which we already know. For instance, for the $\ell=2$ fundamental
mode, we have $r_{\rm s}\,\omega_{2,0}^{\text{Sch}}=0.74734-0.17792\,i$, and
hence

$\displaystyle\omega^{\rm DHOST}_{2,0}=\frac{0.74734-0.17792\,i}{r_{\rm
s}(1+\zeta)^{3/2}}.$ (46)

Equation (45) shows that, even if we know QNM frequencies for multiple pairs
of $(\ell,n)$ from observations, we can determine only the combination $r_{\rm
s}(1+\zeta)^{3/2}$. In this sense, we conclude that there is a degeneracy
between $r_{\rm s}$ and $\zeta$. Another important consequence is that we can
find the QNM frequencies of the stealth Schwarzschild solutions in the DHOST
theory from those in GR by applying the formula (45). This is consistent with
the result of Mukohyama:2023xyf  where the QNM frequencies have been studied
based on the effective field theory with a timelike scalar profile
Mukohyama:2022enj ; Mukohyama:2022skk  applied to a static and spherically
symmetric black hole background.

Let us now briefly discuss the behavior of the power-law tail in the DHOST
theory, assuming that it exists. It is well known that the power-law tail
dominates the waveform of the black hole perturbations after the damped
oscillation phase in GR Price:1971fb . Now, suppose that the solution to the
wave equation (42) (i.e., the one rewritten in the form of the standard Regge-
Wheeler equation in GR) asymptotically behaves as
$\Psi_{\ell}\sim\tilde{T}^{\mathtt{k}}$ at late time, where $\mathtt{k}$ is a
negative constant. Then, by use of Eq. (40), we find that
$\Psi_{\ell}\sim(1+\zeta)^{-\mathtt{k}/2}\tilde{t}^{\mathtt{k}}\propto\tilde{t}^{\mathtt{k}}$.
Therefore, if the power-law tail exists in the DHOST theory, we expect that
its power would be the same as the one in GR.

Before concluding this subsection, we mention the need for a time-domain
analysis. In GR, when we obtain the time evolution of the black hole
perturbations as a solution of Cauchy problem, the late-time behavior of the
black hole perturbations is dominated by a superposition of the QNMs (and the
power-law tail). On the other hand, in the DHOST (or any other modified
gravity) theories, it is nontrivial whether or not the same thing happens
because the effective metric for perturbations does not coincide with the
background metric in general. Although we neglect matter fields in the present
paper, they exist in reality and their dynamics is determined by the
background metric, provided that they are minimally coupled to gravity.
Therefore, in order to obtain the time evolution of the perturbations, we
should impose initial conditions on a hypersurface which is spacelike with
respect to both the background metric and the effective metric. In GR, for
example, the initial surface is often chosen to be a hypersurface with
constant Killing time. In the present case of DHOST theories, a portion of a
constant-$\tilde{t}$ hypersurface can be timelike with respect to the
background metric for $\zeta>0$. When we impose the initial conditions in the
region where the constant-$\tilde{t}$ hypersurface is spacelike, the late-time
behavior of the perturbations would be dominated by the QMNs with frequencies
$\omega^{\rm{DHOST}}_{\ell,n}$. However, when we impose the initial conditions
in the region where the constant-$\tilde{t}$ hypersurface is timelike, it is
not obvious whether the QNMs dominate the late-time behavior of the
perturbations because $\tilde{t}$ cannot be regarded as a physical time
coordinate in this case. In the next subsection, we show that we can prepare a
hypersurface which is spacelike with respect to both the background metric and
the effective metric in the region where the constant-$\tilde{t}$ hypersurface
is timelike by tilting the constant-$\tilde{t}$ hypersurface in an appropriate
manner.

###  III.5 Initial value problem and excitations of QNMs

As we mentioned in Sec. III.2, for $\zeta>0$, a constant-$\tilde{t}$ surface
is timelike with respect to the background metric in the region $r_{\rm
g}<r<r_{\rm g}^{2}/r_{\rm s}$. Therefore, in order to discuss the time
evolution of the perturbations based on the mater equation (38) in a
physically sensible manner, we need to choose another initial hypersurface
that is spacelike with respect to both the background metric and the effective
metric. Such a formulation of initial value problem has been proposed in
Nakashi:2022wdg , which we adopt in the following. In what follows, we focus
on the case with $\zeta>0$ and study the $\ell=2$ mode for concreteness (and
hence the subscript $\ell$ will be omitted). We analyze the initial value
problem for $\zeta<0$ in Appendix A.

Let us briefly review how we construct an initial surface in the physically
sensible formulation proposed in Nakashi:2022wdg . We introduce new
coordinates so that $\tilde{\mathcal{U}}=a\tilde{u}$ and
$\tilde{\mathcal{V}}=b\tilde{v}$, where $a$ and $b$ are positive constants.

Figure 3:  Schematic picture of the initial surface $\Sigma$ (black dashed
curve) and the surface $\tilde{\Sigma}$ (orange dashed curve) which is
constructed by tilting the constant-$\tilde{t}$ surface. We put an initial
Gaussian wave packet (blue solid curve) on the initial surface $\Sigma$. We
require that (a) the initial surface $\Sigma$ coincides with the surface
$\tilde{\Sigma}$ in the region $S$ (red solid curve) within the numerical
domain $D$ (green shaded region). We also require that (b) the initial
conditions have a compact support in the region $S\cap D$, and hence the field
vanishes outside the numerical domain (gray shaded region). We further assume
that the derivative of the field in the direction perpendicular to $\Sigma$ is
zero on the initial surface.

Figure 3 shows a schematic picture of the initial surface and the numerical
domain. The left panel shows our numerical setup in the Penrose diagram of the
Schwarzschild spacetime, while the right panel shows it in a diagram in which
the characteristic curves of the odd-parity perturbations are depicted by 45-
and 135-degree straight lines. By adjusting the constants $a$ and $b$, we can
make a hypersurface of constant $\tilde{\mathcal{U}}+\tilde{\mathcal{V}}$ be
spacelike in the region $r>r_{\rm B}$ for some $r_{\rm B}<r_{\rm g}^{2}/r_{\rm
s}$. We call the constant-$(\mathcal{\tilde{U}}+\mathcal{\tilde{V}})$ surface
$\tilde{\Sigma}$. Let $S$ be the region where the hypersurface
$\tilde{\Sigma}$ is spacelike and let $\Sigma$ denote a spacelike hypersurface
on which we impose initial conditions and the numerical domain $D$. We impose
the following requirements on the hypersurface $\Sigma$ and the initial
conditions:

  1. (a)

The initial surface $\Sigma$ coincides with $\tilde{\Sigma}$ in the region
$S\cap D$.

  2. (b)

The initial conditions have a compact support in the region $S\cap D$.

Since the numerical domain is a part of the causal future of the region $S$
determined by the characteristic curves of the odd-parity perturbations, in
the numerical domain, imposing initial conditions on $\Sigma$ corresponds to
imposing initial conditions on $\tilde{\Sigma}$ under the requirement (a).
Also, under the requirement (a), we can regard
$\tilde{\mathcal{U}}+\tilde{\mathcal{V}}$ as a physical time in the numerical
domain. The requirement (b) allows us to obtain the time evolution as follows.
First, we can obtain the time evolution in the region I in the right panel of
Fig. 3 from the initial data given in the region $S\cap D$. Then, when we
study the time evolution in the regions II and III, we can use the requirement
(b) to set $\Psi=0$ on both the right boundary of region II and the left
boundary of region III. Once we obtain the solution in the regions II and III,
it is straightforward to compute the time evolution in the region IV. Thus, we
can obtain the time evolution of the field in the whole numerical domain.

We consider a Gaussian wave packet as the initial field profile:

$\displaystyle\Psi|_{\Sigma}=\Psi(C_{0}-\tilde{\mathcal{V}},\tilde{\mathcal{V}})|_{\tilde{\Sigma}}=e^{-\frac{1}{2}\left(\frac{\tilde{\mathcal{V}}-\tilde{\mathcal{V}}_{0}}{\sigma}\right)^{2}},$
(47)

where $\sigma$ and $\tilde{\mathcal{V}}_{0}$ are the width of the Gaussian
wave packet and its peak position, respectively. It should be noted that we
truncate the Gaussian profile in a finite region in our actual computations so
that the initial data have a compact support within the region $S\cap D$. We
recall that $\tilde{\mathcal{U}}+\tilde{\mathcal{V}}$ takes a constant value
(which we denote by $C_{0}$) on the initial surface within the numerical
domain thanks to the requirement (a), and hence it makes sense to define the
initial data as in Eq. (47). We choose $\sigma$ and $\tilde{\mathcal{V}}_{0}$
so that the support of the initial field profile overlaps with the region
$r_{\rm B}<r<r_{\rm g}^{2}/r_{\rm s}$, where the surface of constant
$\tilde{\mathcal{U}}+\tilde{\mathcal{V}}$ is spacelike and the
constant-$\tilde{t}$ surface is timelike (see Fig. 4).

Figure 4:  Schematic picture of the initial field profiles. The cyan and the
orange curves are the Gaussian wave packets with $\sigma=r_{\rm s}$ and
$\sigma=0.1r_{\rm s}$, respectively. The horizontal axis is the value of
$(\tilde{\mathcal{U}}-\tilde{\mathcal{V}})/r_{\rm s}$ on the initial surface
$\Sigma$. In practice, we truncate the Gaussian function to have a compact
support within the region $S\cap D$ and choose the center of the wave packet
so that its support overlaps with the region $r_{\rm B}<r<r_{\rm g}^{2}/r_{\rm
s}$. The black curve is the schematic plot of the effective potential
$V(\tilde{\mathcal{U}},\tilde{\mathcal{V}})$ on the initial surface, which is
meant to show the position of the potential peak (and hence its height does
not have any particular meaning).

Also, regarding the initial condition for the derivative, we impose
$(\partial_{\tilde{\mathcal{U}}}+\partial_{\tilde{\mathcal{V}}})\Psi|_{\Sigma}=0$.

Let us now explain how we solve the master equation (38) under the initial
conditions mentioned above. Expressing the master equation (38) in terms of
$\tilde{\mathcal{U}}$ and $\tilde{\mathcal{V}}$, we have

$\displaystyle-4\frac{\partial^{2}\Psi}{\partial\tilde{\mathcal{U}}\partial\tilde{\mathcal{V}}}=\frac{V(\tilde{\mathcal{U}},\tilde{\mathcal{V}})}{ab}\Psi.$
(48)

We discretize the coordinates $\tilde{\mathcal{U}}$ and $\tilde{\mathcal{V}}$
as $\\{\tilde{\mathcal{U}}_{i},\tilde{\mathcal{V}}_{j}\\}$ where
$i,j=0,1,2,\cdots$. Note that the grid width $h$ is assumed to be uniform:
$h=\tilde{\mathcal{U}}_{i+1}-\tilde{\mathcal{U}}_{i}=\tilde{\mathcal{V}}_{j+1}-\tilde{\mathcal{V}}_{j}$.
Also, we introduce shorthand notations
$\Psi_{i,j}=\Psi(\tilde{\mathcal{U}}_{i},\tilde{\mathcal{V}}_{j})$ and
$V_{i,j}=V(\tilde{\mathcal{U}}_{i},\tilde{\mathcal{V}}_{j})$. Then, we apply
the discretization scheme introduced by Gundlach:1993tp , and the master
equation (48) is simply discretized as

$\displaystyle\Psi_{i+1,j+1}=\Psi_{i+1,j}+\Psi_{i,j+1}-\Psi_{i,j}-\frac{h^{2}}{8}\frac{V_{i,j}}{ab}\left[\Psi_{i+1,j}+\Psi_{i,j+1}\right]+\mathcal{O}(h^{4}).$
(49)

The initial field profile (47) can be implemented as

$\displaystyle\Psi_{i,j}=e^{-\frac{1}{2}\left(\frac{\tilde{\mathcal{V}}_{j}-\tilde{\mathcal{V}}_{0}}{\sigma}\right)^{2}}\Theta(\tilde{\mathcal{V}}-\tilde{\mathcal{V}}_{j_{1}})\,\Theta(\tilde{\mathcal{V}}_{j_{2}}-\tilde{\mathcal{V}}),\qquad(\tilde{\mathcal{U}}_{i},\tilde{\mathcal{V}}_{j})\in
S\cap D.$ (50)

Here, to make the truncation explicit, we have inserted the step functions
(denoted by $\Theta$) so that $\Psi$ is nonvanishing only for
$\tilde{\mathcal{V}}_{j_{1}}\leq\tilde{\mathcal{V}}\leq\tilde{\mathcal{V}}_{j_{2}}$
on the initial surface. Also, the condition
$(\partial_{\tilde{\mathcal{U}}}+\partial_{\tilde{\mathcal{V}}})\Psi|_{\Sigma}=0$
yields

$\displaystyle\Psi_{i,j}=\Psi_{i+1,j+1}+\mathcal{O}(h^{4}),\qquad(\tilde{\mathcal{U}}_{i},\tilde{\mathcal{V}}_{j})\in
S\cap D,$ (51)

and hence we have

$\displaystyle\Psi_{i+1,j+1}=\frac{1}{2}\left[\Psi_{i+1,j}+\Psi_{i,j+1}\right]-\frac{h^{2}}{16}\frac{V_{i,j}}{ab}\left[\Psi_{i+1,j}+\Psi_{i,j+1}\right]+\mathcal{O}(h^{4}),\qquad(\tilde{\mathcal{U}}_{i},\tilde{\mathcal{V}}_{j})\in
S\cap D.$ (52)

Combining the discretized equation (49) as well as the initial conditions (50)
and (52), we can obtain the solution for $\Psi$ in the whole numerical domain.

A caveat should be added here. One may think that a constant-$\phi$ surface
would be a good candidate for the initial surface as it is spacelike with
respect to both the background metric and the effective metric. However, if we
choose a constant-$\phi$ surface as the initial surface, it is nontrivial how
to implement initial conditions in our numerical scheme which is based on a
double-null grid, since a constant-$\phi$ surface cannot be described by a
linear function of the null coordinates. This is the reason why we have chosen
the initial surface $\Sigma$ as above.*3*3*3If one would like to choose a
constant-$\phi$ surface as the initial surface, then one needs a numerical
scheme that is more suitable for solving the differential equation based on
the constant-$\phi$ foliation.

Figure 5: The time evolution of the odd-parity perturbations for $\zeta=0.6$
(green curve, top), $\zeta=0.05$ (red curve, middle), and $\zeta=0$ (blue
curve, bottom). The width of the initial Gaussian wave packet is
$\sigma=0.1r_{\rm s}$.

Figure 5 shows the time evolution of the odd-parity perturbations for
$\zeta=0.6$ (green curve, top), $\zeta=0.05$ (red curve, middle), and
$\zeta=0$, i.e., GR (blue curve, bottom) for the initial Gaussian wave packet
with $\sigma=0.1r_{\rm s}$. The observer is located at $\tilde{x}=40\>r_{\rm
s}$. The initial Gaussian wave packet first reaches the observer almost
unscattered, and then the ringdown phase follows, as can be seen in Fig. 5.

In order to confirm that the frequencies in the ringdown phase are QNM
frequencies, we fit the numerical waveform with a superposition of the QNMs.
We introduce the following fitting model $\psi_{N}(\tilde{t})$:

$\displaystyle\psi_{N}(\tilde{t})=\sum_{n=0}^{N}\alpha_{n}e^{-i\left[\mu\,\omega^{\text{Sch}}_{n}(\tilde{t}-\tilde{t}_{\rm
peak})/r_{\rm
s}+\beta_{n}\right]}+c.c.,\qquad\tilde{t}\in[\tilde{t}_{0},\tilde{t}_{\rm
end}],$ (53)

where $n$ labels the overtones and $N$ is the maximum overtone number used in
the fitting. Here, $\alpha_{n}$ and $\beta_{n}$ are real parameters
corresponding to the amplitude and the phase, respectively, and $\mu$ is a
real parameter characterizing the deviation from the QNM frequencies in GR.
Also, $t_{\rm peak}$ denotes the time at which the numerical waveform
$\Psi(\tilde{t})$ takes the maximum value after the initial Gaussian wave
packet passes through the observer. For the fitting analysis, we use the
numerical waveform in the interval $[\tilde{t}_{0},\tilde{t}_{\rm end}]$ where
$\tilde{t}_{0}$ and $\tilde{t}_{\rm end}$ are free parameters satisfying
$\tilde{t}_{\rm peak}\leq\tilde{t}_{0}<\tilde{t}_{\rm end}$. In our fits, we
use the Mathematica function $\mathtt{NonlinearModelFit}$. The amplitude
$\alpha_{n}$ and the phase $\beta_{n}$ are fitting parameters, and we find
best-fit values of these parameters. Note that the parameter $\mu$ is fixed in
the fitting analysis for this section. Once we obtain a best-fit function
$\psi_{N}(\tilde{t})$, we evaluate the goodness of the fit by calculating the
mismatch $\mathcal{M}$ defined by

$\displaystyle\mathcal{M}=1-\frac{\langle\Psi|\psi_{N}\rangle}{\sqrt{\langle\Psi|\Psi\rangle\langle\psi_{N}|\psi_{N}\rangle}},$
(54)

where the scalar product is defined as

$\displaystyle\langle f|g\rangle=\int_{\tilde{t}_{0}}^{\tilde{t}_{\rm
end}}f(\tilde{t})g^{*}(\tilde{t})\;{\rm d}\tilde{t},$ (55)

for arbitrary two complex functions $f$ and $g$, with an asterisk denoting the
complex conjugation. Note that the fitting model $\psi_{N}(\tilde{t})$ and the
numerical waveform $\Psi(\tilde{t})$ are real because we impose real initial
conditions, and hence the complex conjugate in Eq. (55) is of no particular
significance. We have kept the complex conjugate just to follow the convention
in the literature. For the value of $\mu$ which is fixed, we consider the
following two cases:

  1. (i)

$\mu=r_{\rm s}(1+\zeta)^{-3/2}$: We assume the value of $\mu$ as $\mu=r_{\rm
s}\,(1+\zeta)^{-3/2}$ with fixed $\zeta$, and find the best-fit parameters
$\alpha_{n}$ and $\beta_{n}$. This case corresponds to the situation where we
fit the numerical waveform with a superposition of QNMs with the frequencies
$\omega^{\rm DHOST}_{n}$ in the DHOST theory .

  2. (ii)

$\mu=r_{\rm s}$: We assume the value of $\mu$ as $\mu=r_{\rm s}$, and find the
best-fit parameters $\alpha_{n}$ and $\beta_{n}$. This case corresponds to the
situation where we fit the numerical waveform with a superposition of QNMs
with frequencies $\omega^{\text{Sch}}_{n}$ in GR.

Figure 6: The mismatch for $\zeta=0.6$ between the numerical waveform
$\Psi(\tilde{t})$ and the fitting model $\psi_{N}(\tilde{t})$ defined by Eq.
(53). The left panels are calculated with (i) $\mu=r_{\rm s}(1+\zeta)^{-3/2}$,
while the right panels are calculated with (ii) $\mu=r_{\rm s}$. In addition,
the upper panels are the mismatch for $\sigma=r_{\rm s}$, i.e., the wider
initial Gaussian wave packet, while the lower panels are the mismatch for
$\sigma=0.1r_{\rm s}$, i.e., the narrower initial Gaussian wave packet.

Figure 6 shows the mismatch for $\zeta=0.6$. The left panels are the mismatch
calculated with (i) $\mu=r_{\rm s}(1+\zeta)^{-3/2}$ for $\sigma=r_{\rm s}$
(upper panel) and $\sigma=0.1r_{\rm s}$ (lower panel), respectively. As the
number of $N$ increases, the minimum of mismatch decreases. When we fit the
numerical waveform with only the fundamental mode, i.e., $N=0$ (blue curve),
the mismatch gets smaller as $\tilde{t}_{0}$ increases. This implies that the
waveform near the peak time $\tilde{t}_{\rm peak}$ is dominated by the
overtones. Indeed, when we take into account the higher overtones, the value
of $\tilde{t}_{0}$ that minimizes the mismatch gets closer to $\tilde{t}_{\rm
peak}$. The right panels in Fig. 6 show the mismatch calculated with (ii)
$\mu=r_{\rm s}$. Unlike the DHOST fitting (i), for all $N$, the mismatch takes
almost constant values. In particular, we see that the mismatch for $N=0$ does
not decrease at late time in the GR fitting (ii). This reflects the
inconsistency between the numerical waveform and the fitting model: We are now
fitting the waveform for the DHOST theory by a superposition of the QNMs in
GR. Thus, the right panels in Fig. 6 explicitly show that the the QNMs with
the frequencies in GR do not well describe the numerical waveform for the
DHOST theory.

Figure 7: The mismatch for $\zeta=0.05$ between the numerical waveform
$\Psi(\tilde{t})$ and the fitting model $\psi_{N}(\tilde{t})$ defined by Eq.
(53). The values of the parameters in each panel are the same as those in Fig.
6.

Figure 7 shows the mismatch for $\zeta=0.05$. As in the case of $\zeta=0.6$,
when we calculate the mismatch with (i) $\mu=r_{\rm s}(1+\zeta)^{-3/2}$, as
$N$ increases, the minimum of mismatch decreases and the value of
$\tilde{t}_{0}$ at the minimum gets closer to $\tilde{t}_{\rm peak}$. On the
other hand, when we calculate the mismatch with (ii) $r=r_{\rm s}$, it can be
seen that the numerical waveform is not well described with the QNMs in GR.

From these results, we conclude that the superposition of the QNMs in the
DHOST theory (45) is consistent with the numerical waveform and the QNMs are
excited in the physically sensible initial value problem. We have also
performed the fitting analysis keeping $\mu$ unfixed and found that the best-
fit value of $\mu$ is consistent with that of the DHOST theory, i.e.,
$\mu=r_{\rm s}(1+\zeta)^{-3/2}$.

##  IV Stealth Schwarzschild-de Sitter solutions

In this section, we consider the stealth Schwarzschild-dS profile as the
background solution. The background metric is given by the Schwarzschild-dS
metric:

$\displaystyle
A(r)=B(r)=-\frac{\Lambda}{3r}\left(r^{3}-\frac{3}{\Lambda}r+\frac{3r_{\rm
s}}{\Lambda}\right)\eqqcolon-\frac{\Lambda}{3r}\Delta(r),$ (56)

with $r_{\rm s}$ and $\Lambda$ being positive constants. Since $\Delta(r)$ is
a cubic polynomial in $r$, it can be factorized as
$\Delta(r)=(r-r_{-})(r-r_{\rm e})(r-r_{\rm c})$. Here, the three roots are
given by

$\displaystyle r_{-}$
$\displaystyle=\frac{2}{\sqrt{\Lambda}}\cos\left[\frac{1}{3}\cos^{-1}\left(-\frac{3r_{\rm
s}\sqrt{\Lambda}}{2}\right)+\frac{2\pi}{3}\right],$ (57) $\displaystyle r_{\rm
e}$
$\displaystyle=\frac{2}{\sqrt{\Lambda}}\cos\left[\frac{1}{3}\cos^{-1}\left(-\frac{3r_{\rm
s}\sqrt{\Lambda}}{2}\right)+\frac{4\pi}{3}\right],$ (58) $\displaystyle r_{\rm
c}$
$\displaystyle=\frac{2}{\sqrt{\Lambda}}\cos\left[\frac{1}{3}\cos^{-1}\left(-\frac{3r_{\rm
s}\sqrt{\Lambda}}{2}\right)\right],$ (59)

respectively. We note that the three roots are real if

$\displaystyle r_{\rm s}\sqrt{\Lambda}<\frac{2}{3},$ (60)

is satisfied, and we have labeled these roots so that $r_{-}<0<r_{\rm
e}<r_{\rm c}$. Therefore, the event horizon is located at $r=r_{\rm e}$ and
the cosmological horizon is located at $r_{\rm c}$. For small $r_{\rm
s}\sqrt{\Lambda}$, the three roots above are expanded as

$\displaystyle r_{-}$
$\displaystyle=-\sqrt{\frac{3}{\Lambda}}\left[1+\frac{r_{\rm
s}}{2}\sqrt{\frac{\Lambda}{3}}-\frac{r_{\rm s}^{2}\Lambda}{8}+\frac{r_{\rm
s}^{3}}{6}\sqrt{\frac{\Lambda^{3}}{3}}+{\cal O}(r_{\rm
s}^{4}\Lambda^{2})\right],$ (61) $\displaystyle r_{\rm e}$
$\displaystyle=r_{\rm s}\left[1+\frac{1}{3}r_{\rm s}^{2}\Lambda+{\cal
O}(r_{\rm s}^{4}\Lambda^{2})\right],$ (62) $\displaystyle r_{\rm c}$
$\displaystyle=\sqrt{\frac{3}{\Lambda}}\left[1-\frac{r_{\rm
s}}{2}\sqrt{\frac{\Lambda}{3}}-\frac{r_{\rm s}^{2}\Lambda}{8}-\frac{r_{\rm
s}^{3}}{6}\sqrt{\frac{\Lambda^{3}}{3}}+{\cal O}(r_{\rm
s}^{4}\Lambda^{2})\right].$ (63)

###  IV.1 Effective metric

For the stealth Schwarzschild-dS solutions, we can define the effective metric
in the same manner as in the case of stealth Schwarzschild solutions. From the
general expression (29) for the effective metric $Z_{IJ}$, the
$\tilde{t}\tilde{t}$-component can be read off as

$\displaystyle
Z_{\tilde{t}\tilde{t}}=\frac{1}{2F_{2}(1+\zeta)^{2}r^{2}}\left[1-\frac{(1+\zeta)r_{\rm
s}}{r}-\frac{(1+\zeta)\Lambda}{3}r^{2}\right].$ (64)

Introducing a parameter $\Lambda_{\rm g}=(1+\zeta)\Lambda$, we have

$\displaystyle Z_{\tilde{t}\tilde{t}}=-\frac{\Lambda_{\rm
g}}{6F_{2}(1+\zeta)^{2}r^{3}}\left(r^{3}-\frac{3}{\Lambda_{\rm
g}}r+\frac{3r_{\rm g}}{\Lambda_{\rm g}}\right)\eqqcolon-\frac{\Lambda_{\rm
g}}{6F_{2}(1+\zeta)^{2}r^{3}}\Delta_{\rm g}(r).$ (65)

Therefore, for the odd-parity perturbations, the locations of the Killing
horizons are determined by the roots for $\Delta_{\rm g}(r)=0$. The roots are
given by

$\displaystyle\tilde{r}_{-}$ $\displaystyle=\frac{2}{\sqrt{\Lambda_{\rm
g}}}\cos\left[\frac{1}{3}\cos^{-1}\left(-\frac{3r_{\rm g}\sqrt{\Lambda_{\rm
g}}}{2}\right)+\frac{2\pi}{3}\right],$ (66) $\displaystyle\tilde{r}_{\rm e}$
$\displaystyle=\frac{2}{\sqrt{\Lambda_{\rm
g}}}\cos\left[\frac{1}{3}\cos^{-1}\left(-\frac{3r_{\rm g}\sqrt{\Lambda_{\rm
g}}}{2}\right)+\frac{4\pi}{3}\right],$ (67) $\displaystyle\tilde{r}_{\rm c}$
$\displaystyle=\frac{2}{\sqrt{\Lambda_{\rm
g}}}\cos\left[\frac{1}{3}\cos^{-1}\left(-\frac{3r_{\rm g}\sqrt{\Lambda_{\rm
g}}}{2}\right)\right].$ (68)

We note that these three roots are real if $r_{\rm g}\sqrt{\Lambda_{\rm
g}}<2/3$, i.e.,

$\displaystyle r_{\rm s}\sqrt{\Lambda}<\frac{2}{3(1+\zeta)^{3/2}},$ (69)

is satisfied, and we have labeled these roots so that
$\tilde{r}_{-}<0<\tilde{r}_{\rm e}<\tilde{r}_{\rm c}$.

In the present paper, we consider the spacetime region satisfying $\Delta_{\rm
g}(r)>0$, i.e., the static region for the odd-parity perturbations. The static
region associated with the background metric is defined by $\Delta(r)>0$, and
the relation to the static region for the effective metric depends on the
value of $\zeta$ and $\Lambda$. For $\zeta>0$ and $r_{\rm
s}\sqrt{\Lambda}<2(1+\zeta)^{-3/2}/3$, both the effective metric and the
background metric have the static region. However, for
$2(1+\zeta)^{-3/2}/3<r_{\rm s}\sqrt{\Lambda}<2/3$, only the background metric
has the static region.

Figure 8:  Relation between the static region associated with the effective
metric and the one associated with the background metric for $\zeta>0$. The
former (green shaded region) exists for $r_{\rm
s}\sqrt{\Lambda}<2(1+\zeta)^{-3/2}/3$, while the latter (yellow shaded region)
exists for $r_{\rm s}\sqrt{\Lambda}<2/3$.

Figure 8 shows the schematic picture of the relation between the static
regions for the effective metric and the background metric in the $\zeta>0$
case. The solid blue and orange curves respectively correspond to the event
horizon and the cosmological horizon of the background spacetime, while the
dashed blue and orange curves correspond to the Killing horizons for the odd-
parity perturbations. The green shaded region is the static region for both
the effective and background metrics, while the yellow shaded region is the
static region for the background metric only.

###  IV.2 Characters of a constant-$\tilde{t}$ surface

We examine the characters of a constant-$\tilde{t}$ surface. For the stealth
Schwarzschild-dS solutions, the coordinate $\tilde{t}$ defined by (19) reads

$\displaystyle\tilde{t}=t-\int\frac{\zeta(3r)^{3/2}\sqrt{3r_{\rm s}+\Lambda
r^{3}}}{(1+\zeta)\Lambda^{2}\Delta(r)\Delta_{\rm g}(r)}\,{\rm d}r,$ (70)

though we cannot express it in a simple analytic form unlike the case of the
stealth Schwarzschild solutions. We investigate the global structure of a
constant-$\tilde{t}$ surface in the background spacetime. To this end, we
consider a vector field $\partial_{\mu}\tilde{t}$ which is normal to the
constant-$\tilde{t}$ surface. For the background Schwarzschild-dS metric, the
norm of $\partial_{\mu}\tilde{t}$ is given by

$\displaystyle\bar{g}^{\mu\nu}\partial_{\mu}\tilde{t}\,\partial_{\nu}\tilde{t}=\frac{3r}{\Lambda\Delta_{\rm
g}(r)}\left(r^{3}-\frac{3}{(1+\zeta)^{2}\Lambda}r+\frac{3r_{\rm
s}}{\Lambda}\right).$ (71)

Since we focus only on the region in which $\Delta_{\rm g}(r)>0$, the sign of
$\bar{g}^{\mu\nu}\partial_{\mu}\tilde{t}\,\partial_{\nu}\tilde{t}$ is
determined by the sign of the function in the parentheses in Eq. (71).

Figure 9:  Typical plots of constant-$\tilde{t}$ surfaces in the Penrose
diagram of the Schwarzschild-dS spacetime. The black curves represent the
constant-$\tilde{t}$ surfaces. The constant-$\tilde{t}$ surfaces are spacelike
in the yellow shaded region. For (A-1) $\zeta>0$, when the parameter $\zeta$
satisfies Eq. (72), the constant-$\tilde{t}$ surface can be spacelike in a
finite region. For (A-2) $\zeta>0$ when the parameter $\zeta$ violates Eq.
(72), the constant-$\tilde{t}$ surface is always timelike. For (B) $\zeta<0$,
the constant-$\tilde{t}$ surface is always spacelike.

In Fig. 9, we show typical plots of constant-$\tilde{t}$ surfaces in the
Penrose diagram of the Schwarzschild-dS spacetime.*4*4*4 In Fig. 9 and Fig.
10, in the $\zeta<0$ case, a constant-$\tilde{t}$ surface and a characteristic
curves apparently become null near both the event horizon and the cosmological
horizon, but in fact these are spacelike everywhere. This behavior is an
artifact caused by the coordinate system we have used. For the stealth
Schwarzschild-dS solutions, we have defined different double null coordinates
in each block of the Penrose diagram and drawn the curves for each block.
Then, we have glued the diagrams of each block using the method proposed in
Walker1970 . On the other hand, for the stealth Schwarzschild solutions, we
can define a single coordinate system which covers the whole spacetime in a
simple analytic form, and hence the coordinate $\tilde{t}$ and the
characteristic curves are written in an analytic form, which has led to the
smooth curves in Fig. 1 and Fig. 2.  The constant-$\tilde{t}$ surfaces are
spacelike in the yellow shaded region. For $\zeta>0$, the constant-$\tilde{t}$
surface has a spacelike region within $\tilde{r}_{\rm e}<r<\tilde{r}_{\rm c}$
if

$\displaystyle r_{\rm s}\sqrt{\Lambda}<\frac{2}{3(1+\zeta)^{3}},$ (72)

is satisfied [see (A-1) in Fig. 9]. If the condition (72) is violated for
$\zeta>0$, the constant-$\tilde{t}$ surface is always timelike [see (A-2) in
Fig. 9]. On the other hand, for $\zeta<0$, the constant-$\tilde{t}$ surface is
always spacelike [see (B) in Fig. 9].

###  IV.3 Characteristic curves

As in the case of the stealth Schwarzschild solutions, the odd-parity
perturbations propagate along a constant-$\tilde{u}$ curve or a
constant-$\tilde{v}$ curve. Here, we investigate the vector fields which are
normal to the constant-$\tilde{u}$ curve and the constant-$\tilde{v}$ curve:
$\partial_{\mu}\tilde{u}$ and $\partial_{\mu}\tilde{v}$. For the background
Schwarzschild-dS metric, the norm of these vector fields are given by

$\displaystyle\bar{g}^{\mu\nu}\partial_{\mu}\tilde{u}\,\partial_{\nu}\tilde{u}$
$\displaystyle=\frac{3\zeta r}{\Lambda_{\rm g}^{2}\Delta_{\rm
g}(r)^{2}}\left(\sqrt{3r_{\rm g}+\Lambda_{\rm g}r^{3}}+\sqrt{3r}\right)^{2},$
(73)
$\displaystyle\bar{g}^{\mu\nu}\partial_{\mu}\tilde{v}\,\partial_{\nu}\tilde{v}$
$\displaystyle=\frac{3\zeta r}{\Lambda_{\rm g}^{2}\Delta_{\rm
g}(r)^{2}}\left(\sqrt{3r_{\rm g}+\Lambda_{\rm g}r^{3}}-\sqrt{3r}\right)^{2}.$
(74)

The sign of each norm is determined by the sign of the parameter $\zeta$. The
characteristic curves are timelike for $\zeta>0$, while they are spacelike for
$\zeta<0$. That is, for $\zeta<0$, the odd-parity perturbations become
superluminal. Figure 10 shows the characteristic curves of the odd-parity
perturbations in the Penrose diagram of the Schwarzschild-dS spacetime.

Figure 10: The characteristic curves for (A) $\zeta>0$ and (B) $\zeta<0$ in
the Penrose diagram of the Schwarzschild-dS spacetime. The red curves and the
blues curves represent the constant-$\tilde{v}$ curves and the
constant-$\tilde{u}$ curves, respectively.

###  IV.4 Equation of motion and QNM frequencies

In order to express the equation of motion in the form of a two-dimensional
wave equation, we introduce the generalized tortoise coordinate:

$\displaystyle\tilde{x}=\frac{-3}{\Lambda\sqrt{1+\zeta}}\left(\frac{\tilde{r}_{-}\ln\left|r-\tilde{r}_{-}\right|}{(\tilde{r}_{\rm
c}-\tilde{r}_{-})(\tilde{r}_{\rm e}-\tilde{r}_{-})}+\frac{\tilde{r}_{\rm
e}\ln\left|r-\tilde{r}_{\rm e}\right|}{(\tilde{r}_{\rm e}-\tilde{r}_{\rm
c})(\tilde{r}_{\rm e}-\tilde{r}_{-})}+\frac{\tilde{r}_{\rm
c}\ln\left|r-\tilde{r}_{\rm c}\right|}{(\tilde{r}_{\rm c}-\tilde{r}_{\rm
e})(\tilde{r}_{\rm c}-\tilde{r}_{-})}\right),$ (75)

up to an integration constant. We note that $\tilde{x}\to-\infty$ as
$r\to\tilde{r}_{\rm e}$ and $\tilde{x}\to\infty$ as $r\to\tilde{r}_{\rm c}$.
In terms of $\tilde{t}$ and $\tilde{x}$, the equation of motion is written as
follows:

$\displaystyle\left[\frac{\partial^{2}}{\partial\tilde{x}^{2}}-\frac{\partial^{2}}{\partial\tilde{t}^{2}}-V_{\ell}(\tilde{x})\right]\Psi_{\ell}=0,$
(76)

with

$\displaystyle V_{\ell}(\tilde{x})=\frac{1}{1+\zeta}\left(1-\frac{r_{\rm
g}}{r}-\frac{\Lambda_{\rm
g}}{3}r^{2}\right)\left[\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{\rm
g}}{r^{3}}\right].$ (77)

As in the case of the stealth Schwarzschild solutions, we introduce the
rescaled coordinates as follows:

$\displaystyle\tilde{T}=\frac{\tilde{t}}{\sqrt{1+\zeta}},$ (78)
$\displaystyle\tilde{X}=\frac{\tilde{x}}{\sqrt{1+\zeta}}.$ (79)

With these coordinates, the equation of motion can be rewritten as

$\displaystyle\left[\frac{\partial^{2}}{\partial\tilde{X}^{2}}-\frac{\partial^{2}}{\partial\tilde{T}^{2}}-\tilde{V}_{\ell}(\tilde{X})\right]\Psi_{\ell}=0,$
(80)

with

$\displaystyle\tilde{V}_{\ell}(\tilde{X})=\left(1-\frac{r_{\rm
g}}{r}-\frac{\Lambda_{\rm
g}}{3}r^{2}\right)\left[\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{\rm
g}}{r^{3}}\right].$ (81)

Therefore, even in the case of stealth Schwarzschild-dS solutions, we can
recast the equation of motion into the standard Regge-Wheeler equation
parametrized by $r_{\rm g}$ and $\Lambda_{\rm g}$. For the stealth
Schwarzschild-dS solutions, the QNMs correspond to the modes that are purely
ingoing at $r=\tilde{r}_{\rm e}$ and purely outgoing at $r=\tilde{r}_{\rm c}$.
To determine the QNM frequencies, we consider the ansatz
$\Psi_{\ell}=\psi_{\ell}(\tilde{X})e^{-i\tilde{W}\tilde{T}}$. Let
$\tilde{W}_{\ell,n}(r_{\rm g},\Lambda_{\rm g})$ be the QNM frequencies
obtained by solving Eq. (80) with the ansatz for $\Psi_{\ell}$ and let
$\omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm s},\Lambda)$ be the QNM frequencies
calculated from the standard Regge-Wheeler equation parametrized by $r_{\rm
s}$ and $\Lambda$. Then, the QNM frequencies $\tilde{W}_{\ell,n}(r_{\rm
g},\Lambda_{\rm g})$ can be expressed as

$\displaystyle\tilde{W}_{\ell,n}(r_{\rm g},\Lambda_{\rm g})=\omega^{\text{Sch-
dS}}_{\ell,n}(r_{\rm g},\Lambda_{\rm g}),$ (82)

where $\omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm g},\Lambda_{\rm
g})\coloneqq\omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm s}\to r_{\rm
g},\Lambda\to\Lambda_{\rm g})$. Finally, from the relation between $\tilde{T}$
and $\tilde{t}$, the QNM frequencies of the odd-parity perturbations about the
stealth Schwarzschild-dS solutions in the DHOST theory are given by

$\displaystyle\omega^{\rm DHOST}_{\ell,n}=\frac{\omega^{\text{Sch-
dS}}_{\ell,n}(r_{\rm g},\Lambda_{\rm g})}{\sqrt{1+\zeta}}.$ (83)

That is, once we know the QNM frequencies in the Schwarzschild-dS spacetime
parametrized by $r_{\rm g}$ and $\Lambda_{\rm g}$ in GR, we can find the QNM
frequencies in the DHOST theory by applying the scaling law Eq. (83).

It is worth mentioning that there is a degeneracy between $r_{\rm s}$ and
$\zeta$ as in the case of the stealth Schwarzschild solutions. As we discuss
in Appendix B, the dimensionless QNM frequencies of the Schwarzschild-dS
spacetime in GR, $\Omega^{\text{Sch-dS}}_{\ell,n}\coloneqq r_{\rm
s}\,\omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm s},\Lambda)$, depends only on the
combination $r_{\rm s}^{2}\Lambda$, and hence we can write $\Omega^{\text{Sch-
dS}}_{\ell,n}=\Omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm s}^{2}\Lambda)$. Thus,
the QNM frequencies of the Schwarzschild-dS spacetime can be written as

$\displaystyle\omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm
s},\Lambda)=\frac{\Omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm
s}^{2}\Lambda)}{r_{\rm s}}.$ (84)

As a result, Eq. (83) can be rewritten in terms of $\Omega^{\text{Sch-
dS}}_{\ell,n}$ as follows:

$\displaystyle\omega^{\rm DHOST}_{\ell,n}=\frac{\Omega^{\text{Sch-
dS}}_{\ell,n}(r_{\rm g}^{2}\Lambda_{\rm g})}{r_{\rm
g}\sqrt{1+\zeta}}=\frac{\Omega^{\text{Sch-dS}}_{\ell,n}\left([r_{\rm
s}(1+\zeta)^{3/2}]^{2}\Lambda\right)}{r_{\rm s}(1+\zeta)^{3/2}},$ (85)

where we have used $r_{\rm g}=(1+\zeta)r_{\rm s}$ and $\Lambda_{\rm
g}=(1+\zeta)\Lambda$. This means that, if we fix the value of $\Lambda$,
$\omega^{\text{DHOST}}_{\ell,n}$ depends only on the combination $r_{\rm
s}(1+\zeta)^{3/2}$. This shows that there is a degeneracy between $r_{\rm s}$
and $\zeta$ in the QNM frequencies for the stealth Schwarzschild-dS solutions
as in the case of the stealth Schwarzschild solutions.

For $r_{\rm g}^{2}\Lambda_{\rm g}\ll 1$, there is a perturbative formula for
the QNM frequencies Hatsuda:2023geo . For $\ell=2$ fundamental mode, the
formula is given by

$\displaystyle\Omega^{\text{Sch-dS}}_{2,0}(r_{\rm g}^{2}\Lambda_{\rm
g})=0.74734-0.17792\,i+\frac{9w_{1}}{2}r_{\rm g}^{2}\Lambda_{\rm
g}+\frac{81w_{2}}{8}r_{\rm g}^{4}\Lambda_{\rm g}^{2}+\mathcal{O}(r_{\rm
g}^{6}\Lambda_{\rm g}^{3}),$ (86)

where $w_{1}=-0.18649+0.03720\,i$, $w_{2}=-0.04819+0.01428\,i$. Note that this
formula implicitly assumes $r_{\rm s}\neq 0$ since otherwise the left-hand
side cannot be defined. Therefore, the QNM frequency $\omega^{\rm
DHOST}_{2,0}$ becomes

$\displaystyle\begin{split}\omega^{\rm DHOST}_{2,0}=\frac{1}{r_{\rm
s}(1+\zeta)^{3/2}}\bigg{(}0.74734&-0.17792\,i+\frac{9w_{1}}{2}\left[r_{\rm
s}(1+\zeta)^{3/2}\right]^{2}\Lambda\\\ &+\frac{81w_{2}}{8}\left[r_{\rm
s}(1+\zeta)^{3/2}\right]^{4}\Lambda^{2}+\mathcal{O}(r_{\rm
s}^{6}\Lambda^{3})\bigg{)}.\end{split}$ (87)

This explicitly shows that there is the degeneracy between $r_{\rm s}$ and
$\zeta$ when we fix the value of $\Lambda$.

##  V Summary and discussions

We have investigated the odd-parity perturbations about stealth Schwarzschild
solutions and stealth Schwarzschild-de Sitter solutions with a linearly time-
dependent scalar field in a subclass of DHOST theories, for which the
deviation from general relativity is controlled by a single parameter $\zeta$.
We have derived the effective metric for the odd-parity perturbations and
analyzed the characteristic curves. We have also shown that the Killing
horizon(s) of the effective metric differs from that of the background metric.
For $\zeta<0$ case, the odd-parity perturbations can be superluminal and hence
can escape from the region inside the Schwarzschild radius of the background
metric, as demonstrated in Appendix A. We have derived the master equation for
the odd-parity perturbations in the form of a two-dimensional wave equation,
which can be expressed in the form of the standard Regge-Wheeler equation in
GR with the rescaled black hole mass $r_{\rm g}$ (and the rescaled
cosmological constant $\Lambda_{\rm g}$ in the case of the stealth
Schwarzschild-de Sitter solutions). We have computed the QNM frequencies for
both the stealth Schwarzschild solutions and the stealth Schwarzschild-dS
solutions. In both cases, we have found that the QNM frequencies can be given
by a simple scaling of those in GR. In particular, we have shown that there is
a degeneracy between the black hole mass $r_{\rm s}$ and $\zeta$ in the QNM
frequencies.

We have also solved an initial value problem for the odd-parity perturbations
about the stealth Schwarzschild solutions employing the physically sensible
formulation of the initial value problem proposed in Nakashi:2022wdg . We have
defined a spacelike hypersurface $\Sigma$ in the following manner: We have
constructed a hypersurface $\tilde{\Sigma}$ by slightly tilting the
constant-$\tilde{t}$ surface. We have defined the region $S$ where the surface
$\tilde{\Sigma}$ is spacelike. Furthermore, we have required that (a) the
initial surface $\Sigma$ coincides with $\tilde{\Sigma}$ in the region $S$
within the numerical domain, and that (b) the initial conditions have a
compact support in the region $S$ within the numerical domain. We have
analyzed the time evolution of a initial Gaussian wave packet. We have
confirmed that the damped oscillation phase (ringdown phase) appears. We have
found that a superposition of the QNMs in the DHOST theory is consistent with
the numerical waveform through the fitting analysis. In particular, we have
calculated the mismatch between the numerical waveform and the superposition
of the QNMs in the DHOST theory and found that the minimum of the mismatch
decreases and gets closer to the waveform peak when the overtones are taken
into account. On the other hand, we have also confirmed that a superposition
of the QNMs in GR does not well describe the numerical waveform. From these
results, we conclude that the QNMs in the DHOST theory are excited in the
physically sensible initial value problem.

We note that the perturbations about the stealth solution in DHOST theories
would be strongly coupled Babichev:2018uiw ; deRham:2019gha ;
Motohashi:2019ymr ; Takahashi:2021bml . A possible way out of this problem is
to incorporate the scordatura term Motohashi:2019ymr . However, as discussed
in Sec. II.3, we expect that the scordatura term would not lead to a
qualitative change in our results on the odd-parity perturbations. This is
essentially because the strong coupling problem comes from the vanishing sound
speed of the mode corresponding to the scalar degree of freedom, which belongs
to the even-parity sector. Along this line of thought, it would be intriguing
to study the even-parity sector to see how the effect of the scordatura term
shows up. It should also be noted that the effect of modified gravity
completely disappears in the odd-parity sector when $\zeta=0$, and hence the
study of odd-parity perturbations alone cannot tell the difference from
general relativity. This is another motivation to study the even-parity
sector. We hope to come back to this issue in a future publication.

## Acknowledgments

This work was supported by JSPS KAKENHI Grant Nos. JP22K03626 (M.K.),
JP22K03639 (H.M.), JP22KJ1646 (K.T.), and JP23K13101 (K.T.) from the Japan
Society for the Promotion of Science.

##  Appendix A Initial value problem for negative $\zeta$

In Sec. III.5, we have analyzed the initial value problem for the odd-parity
perturbations about the stealth Schwarzschild solutions in $\zeta>0$ case and
shown that the waveform of the odd-parity perturbations is well described by a
superposition of the QNMs in the DHOST theory. In this appendix, we analyze
the case of $\zeta<0$. A remarkable property of the odd-parity perturbation
for $\zeta<0$ is that the perturbation is superluminal in the whole spacetime.
This implies that the odd-parity perturbations can escape from inside the
Schwarzschild radius $r_{\rm s}$. As we mentioned in Sec. III.2, a
constant-$\tilde{t}$ surface is always spacelike, and hence we choose it as
the initial surface, i.e., we set $a=b=1$. Furthermore, we choose $\sigma$ and
$\tilde{\mathcal{V}}_{0}$ so that the initial field profile has its support
inside $r_{\rm s}$.

Figure 11:  The evolution of the odd-parity perturbation in the
$(\tilde{\mathcal{U}},\tilde{\mathcal{V}})$-space (left panel) and the
waveform of the perturbation observed by an observer at $\tilde{x}=40r_{\rm
s}$ (right panel) for $\zeta<0$. In the left panel, the cyan solid curve is
the initial field profile and the black solid line is the location of the
Schwarzschild radius $r_{\rm s}$. The left panel shows that the odd-parity
perturbations can escape from inside the Schwarzschild radius. The right panel
shows that the damped oscillation phase also shows up for $\zeta<0$.

Figure 11 shows the evolution of the odd-parity perturbation for $\zeta=-0.5$
in the $(\tilde{\mathcal{U}},\tilde{\mathcal{V}})$-space (left panel) and the
waveform observed by an observer at $\tilde{x}=40r_{\rm s}$ (right panel). In
the left panel, the cyan solid curve is the initial field profile, which has
the compact support inside $r_{\rm s}$ depicted by the black solid line. The
left panel explicitly shows that the odd-parity perturbation escapes from
inside the Schwarzschild radius $r_{\rm s}$. The right panel shows that the
damped oscillation phase also shows up in the $\zeta<0$ case. We note that the
damped oscillation phase can be well fitted by Eq. (53) with $\mu=r_{\rm
s}(1+\zeta)^{-3/2}$.

##  Appendix B QNM frequencies of the Schwarzschild-dS spacetime in GR

Here, we briefly mention a property of the QNMs of the Schwarzschild-dS
spacetime in GR. The Regge-Wheeler equation for the Schwarzschild-dS solution
in GR can be written as

$\displaystyle\left[\frac{\partial^{2}}{\partial
x^{2}}-\frac{\partial^{2}}{\partial t^{2}}-V_{\ell}(x)\right]\Psi_{\ell}=0,$
(88)

with the effective potential given by

$\displaystyle V_{\ell}(x)=\left(1-\frac{r_{\rm
s}}{r}-\frac{\Lambda}{3}r^{2}\right)\left[\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{\rm
s}}{r^{3}}\right].$ (89)

Substituting the ansatz $\Psi_{\ell}=\psi_{\ell}(x)e^{-i\omega t}$, we have

$\displaystyle\left[\frac{\partial^{2}}{\partial
x^{2}}+\omega^{2}-V_{\ell}(x)\right]\psi_{\ell}(x)=0.$ (90)

In terms of the dimensionless coordinates $\hat{r}\coloneqq r/r_{\rm s}$ and
$\hat{x}\coloneqq x/r_{\rm s}$, the above equation takes the form

$\displaystyle\left[\frac{\partial^{2}}{\partial\hat{x}^{2}}+(r_{\rm
s}\omega)^{2}-r_{\rm s}^{2}V_{\ell}(x)\right]\psi_{\ell}(x)=0.$ (91)

Note that the effective potential written in terms of the dimensionless
coordinates reads

$\displaystyle r_{\rm s}^{2}V_{\ell}(x)=\left(1-\frac{1}{\hat{r}}-\frac{r_{\rm
s}^{2}\Lambda}{3}\hat{r}^{2}\right)\left[\frac{\ell(\ell+1)}{\hat{r}^{2}}-\frac{3}{\hat{r}^{3}}\right],$
(92)

where $r_{\rm s}$ and $\Lambda$ show up only in the combination $r_{\rm
s}^{2}\Lambda$. As a result, the QNM frequencies normalized by $r_{\rm s}$
should depend only on $r_{\rm s}^{2}\Lambda$, which we write
$\Omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm s}^{2}\Lambda)$. Therefore, the
(dimensionful) QNM frequencies $\omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm
s},\Lambda)$ can be expressed as

$\displaystyle\omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm
s},\Lambda)=\frac{\Omega^{\text{Sch-dS}}_{\ell,n}(r_{\rm
s}^{2}\Lambda)}{r_{\rm s}},$ (93)

For $r_{\rm s}^{2}\Lambda\ll 1$, there is a perturbative formula for the QNM
frequencies Hatsuda:2023geo . According to the formula, for $\ell=2$
fundamental mode, the dimensionless QNM frequency $\Omega^{\text{Sch-
dS}}_{2,0}(r_{\rm s}^{2}\Lambda)$ is given by

$\displaystyle\Omega^{\text{Sch-dS}}_{2,0}(r_{\rm
s}^{2}\Lambda)=0.74734-0.17792\,i+\frac{9w_{1}}{2}r_{\rm
s}^{2}\Lambda+\frac{81w_{2}}{8}r_{\rm s}^{4}\Lambda^{2}+\mathcal{O}(r_{\rm
s}^{6}\Lambda^{3}),$ (94)

where $w_{1}=-0.18649+0.03720\,i$, $w_{2}=-0.04819+0.01428\,i$. This is
consistent with the fact that the dimensionless QNM frequencies
$\Omega^{\text{Sch-dS}}_{\ell,n}$ depend only on $r_{\rm s}^{2}\Lambda$.

## References

  * (1) K. Koyama, _Cosmological Tests of Modified Gravity_ , _Rept. Prog. Phys._ 79 (2016) 046902 [1504.04623].
  * (2) P. G. Ferreira, _Cosmological Tests of Gravity_ , _Ann. Rev. Astron. Astrophys._ 57 (2019) 335 [1902.10503].
  * (3) S. Arai et al., _Cosmological gravity probes: Connecting recent theoretical developments to forthcoming observations_ , _PTEP_ 2023 (2023) 072E01 [2212.09094].
  * (4) A. Buonanno, G. B. Cook and F. Pretorius, _Inspiral, merger and ring-down of equal-mass black-hole binaries_ , _Phys. Rev. D_ 75 (2007) 124018 [gr-qc/0610122].
  * (5) C. Brans and R. H. Dicke, _Mach’s principle and a relativistic theory of gravitation_ , _Phys. Rev._ 124 (1961) 925.
  * (6) G. W. Horndeski, _Second-order scalar-tensor field equations in a four-dimensional space_ , _Int. J. Theor. Phys._ 10 (1974) 363.
  * (7) C. Deffayet, X. Gao, D. A. Steer and G. Zahariade, _From k-essence to generalised Galileons_ , _Phys. Rev. D_ 84 (2011) 064039 [1103.3260].
  * (8) T. Kobayashi, M. Yamaguchi and J. Yokoyama, _Generalized G-inflation: Inflation with the most general second-order field equations_ , _Prog. Theor. Phys._ 126 (2011) 511 [1105.5723].
  * (9) R. P. Woodard, _Ostrogradsky’s theorem on Hamiltonian instability_ , _Scholarpedia_ 10 (2015) 32243 [1506.02210].
  * (10) H. Motohashi and T. Suyama, _Third order equations of motion and the Ostrogradsky instability_ , _Phys. Rev. D_ 91 (2015) 085009 [1411.3721].
  * (11) H. Motohashi and T. Suyama, _Quantum Ostrogradsky theorem_ , _JHEP_ 09 (2020) 032 [2001.02483].
  * (12) K. Aoki and H. Motohashi, _Ghost from constraints: a generalization of Ostrogradsky theorem_ , _JCAP_ 08 (2020) 026 [2001.06756].
  * (13) D. Langlois and K. Noui, _Degenerate higher derivative theories beyond Horndeski: evading the Ostrogradski instability_ , _JCAP_ 02 (2016) 034 [1510.06930].
  * (14) H. Motohashi, K. Noui, T. Suyama, M. Yamaguchi and D. Langlois, _Healthy degenerate theories with higher derivatives_ , _JCAP_ 07 (2016) 033 [1603.09355].
  * (15) R. Klein and D. Roest, _Exorcising the Ostrogradsky ghost in coupled systems_ , _JHEP_ 07 (2016) 130 [1604.01719].
  * (16) H. Motohashi, T. Suyama and M. Yamaguchi, _Ghost-free theory with third-order time derivatives_ , _J. Phys. Soc. Jap._ 87 (2018) 063401 [1711.08125].
  * (17) H. Motohashi, T. Suyama and M. Yamaguchi, _Ghost-free theories with arbitrary higher-order time derivatives_ , _JHEP_ 06 (2018) 133 [1804.07990].
  * (18) M. Crisostomi, K. Koyama and G. Tasinato, _Extended Scalar-Tensor Theories of Gravity_ , _JCAP_ 04 (2016) 044 [1602.03119].
  * (19) J. Ben Achour, M. Crisostomi, K. Koyama, D. Langlois, K. Noui and G. Tasinato, _Degenerate higher order scalar-tensor theories beyond Horndeski up to cubic order_ , _JHEP_ 12 (2016) 100 [1608.08135].
  * (20) J. D. Bekenstein, _The Relation between physical and gravitational geometry_ , _Phys. Rev. D_ 48 (1993) 3641 [gr-qc/9211017].
  * (21) J.-P. Bruneton and G. Esposito-Farèse, _Field-theoretical formulations of MOND-like gravity_ , _Phys. Rev. D_ 76 (2007) 124012 [0705.4043].
  * (22) D. Bettoni and S. Liberati, _Disformal invariance of second order scalar-tensor theories: Framing the Horndeski action_ , _Phys. Rev. D_ 88 (2013) 084020 [1306.6724].
  * (23) K. Takahashi, H. Motohashi and M. Minamitsuji, _Invertible disformal transformations with higher derivatives_ , _Phys. Rev. D_ 105 (2022) 024015 [2111.11634].
  * (24) K. Takahashi, _Invertible disformal transformations with arbitrary higher-order derivatives_ , 2307.08814.
  * (25) K. Takahashi, M. Minamitsuji and H. Motohashi, _Generalized disformal Horndeski theories: cosmological perturbations and consistent matter coupling_ , _PTEP_ 2023 (2023) 013E01 [2209.02176].
  * (26) A. Naruko, R. Saito, N. Tanahashi and D. Yamauchi, _Ostrogradsky mode in scalar-tensor theories with higher-order derivative couplings to matter_ , _PTEP_ 2023 (2023) 053E02 [2209.02252].
  * (27) K. Takahashi, R. Kimura and H. Motohashi, _Consistency of matter coupling in modified gravity_ , _Phys. Rev. D_ 107 (2023) 044018 [2212.13391].
  * (28) T. Ikeda, K. Takahashi and T. Kobayashi, _Consistency of higher-derivative couplings to matter fields in scalar-tensor gravity_ , _Phys. Rev. D_ 108 (2023) 044006 [2302.03418].
  * (29) A. De Felice, D. Langlois, S. Mukohyama, K. Noui and A. Wang, _Generalized instantaneous modes in higher-order scalar-tensor theories_ , _Phys. Rev. D_ 98 (2018) 084024 [1803.06241].
  * (30) A. De Felice, S. Mukohyama and K. Takahashi, _Nonlinear definition of the shadowy mode in higher-order scalar-tensor theories_ , _JCAP_ 12 (2021) 020 [2110.03194].
  * (31) A. De Felice, S. Mukohyama and K. Takahashi, _Built-in scordatura in U-DHOST_ , _Phys. Rev. Lett._ 129 (2022) 031103 [2204.02032].
  * (32) K. Takahashi, M. Minamitsuji and H. Motohashi, _Effective description of generalized disformal theories_ , _JCAP_ 07 (2023) 009 [2304.08624].
  * (33) H. Nariai, _On the Green’s Function in an Expanding Universe and Its Role in the Problem of Mach’s Principle_ , _Progress of Theoretical Physics_ 40 (1968) 49.
  * (34) J. O’ Hanlon and B. O. J. Tupper, _Vacuum-field solutions in the Brans-Dicke theory_ , _Nuovo Cim. B_ 7 (1972) 305.
  * (35) J. D. Barrow and K. Maeda, _Extended inflationary universes_ , _Nuclear Physics B_ 341 (1990) 294.
  * (36) C. Romero and A. Barros, _Brans-dicke vacuum solutions and the cosmological constant: A qualitative analysis_ , _General Relativity and Gravitation_ 25 (1993) 491.
  * (37) S. J. Kolitch, _Qualitative analysis of Brans-Dicke universes with a cosmological constant_ , _Annals Phys._ 246 (1996) 121 [gr-qc/9409002].
  * (38) V. B. Johri and K. Desikan, _Cosmological models with constant deceleration parameter in Brans-Dicke theory_ , _Gen. Rel. Grav._ 26 (1994) 1217.
  * (39) S. Giardino, V. Faraoni and A. Giusti, _First-order thermodynamics of scalar-tensor cosmology_ , _JCAP_ 04 (2022) 053 [2202.07393].
  * (40) S. Giardino, A. Giusti and V. Faraoni, _Thermal stability of stealth and de Sitter spacetimes in scalar-tensor gravity_ , _Eur. Phys. J. C_ 83 (2023) 621 [2302.08550].
  * (41) E. Ayón-Beato, C. Martínez and J. Zanelli, _Stealth scalar field overflying a (2+1) black hole_ , _Gen. Rel. Grav._ 38 (2006) 145 [hep-th/0403228].
  * (42) E. Ayón-Beato, C. Martínez, R. Troncoso and J. Zanelli, _Gravitational Cheshire effect: Nonminimally coupled scalar fields may not curve spacetime_ , _Phys. Rev. D_ 71 (2005) 104037 [hep-th/0505086].
  * (43) S. Mukohyama, _Black holes in the ghost condensate_ , _Phys. Rev. D_ 71 (2005) 104019 [hep-th/0502189].
  * (44) D. C. Robinson, _Non-gravitating waves_ , _Gen. Rel. Grav._ 38 (2006) 153.
  * (45) E. Ayón-Beato, M. Hassaïne and M. M. Juárez-Aubry, _Stealths on Anisotropic Holographic Backgrounds_ , 1506.03545.
  * (46) A. Alvarez, C. Campuzano, M. Cruz, E. Rojas and J. Saavedra, _Stealths on $(1+1)$-dimensional dilatonic gravity_, _Gen. Rel. Grav._ 48 (2016) 165 [1611.03022].
  * (47) I. Smolić, _Spacetimes dressed with stealth electromagnetic fields_ , _Phys. Rev. D_ 97 (2018) 084041 [1711.07490].
  * (48) E. Franzin and I. Smolić, _Stationary spacetimes with time-dependent real scalar fields_ , _Class. Quant. Grav._ 38 (2021) 115004 [2101.05816].
  * (49) E. Babichev and C. Charmousis, _Dressing a black hole with a time-dependent Galileon_ , _JHEP_ 08 (2014) 106 [1312.3204].
  * (50) T. Kobayashi and N. Tanahashi, _Exact black hole solutions in shift symmetric scalar–tensor theories_ , _PTEP_ 2014 (2014) 073E02 [1403.4364].
  * (51) E. Babichev and G. Esposito-Farèse, _Cosmological self-tuning and local solutions in generalized Horndeski theories_ , _Phys. Rev. D_ 95 (2017) 024020 [1609.09798].
  * (52) E. Babichev, C. Charmousis, G. Esposito-Farèse and A. Lehébel, _Stability of Black Holes and the Speed of Gravitational Waves within Self-Tuning Cosmological Models_ , _Phys. Rev. Lett._ 120 (2018) 241101 [1712.04398].
  * (53) M. Minamitsuji and H. Motohashi, _Stealth Schwarzschild solution in shift symmetry breaking theories_ , _Phys. Rev. D_ 98 (2018) 084027 [1809.06611].
  * (54) J. Ben Achour and H. Liu, _Hairy Schwarzschild-(A)dS black hole solutions in degenerate higher order scalar-tensor theories beyond shift symmetry_ , _Phys. Rev. D_ 99 (2019) 064042 [1811.05369].
  * (55) H. Motohashi and M. Minamitsuji, _General Relativity solutions in modified gravity_ , _Phys. Lett. B_ 781 (2018) 728 [1804.01731].
  * (56) H. Motohashi and M. Minamitsuji, _Exact black hole solutions in shift-symmetric quadratic degenerate higher-order scalar-tensor theories_ , _Phys. Rev. D_ 99 (2019) 064040 [1901.04658].
  * (57) M. Minamitsuji and J. Edholm, _Black hole solutions in shift-symmetric degenerate higher-order scalar-tensor theories_ , _Phys. Rev. D_ 100 (2019) 044053 [1907.02072].
  * (58) R. C. Bernardo, J. Celestial and I. Vega, _Stealth black holes in shift symmetric kinetic gravity braiding_ , _Phys. Rev. D_ 101 (2020) 024036 [1911.01847].
  * (59) C. Charmousis, M. Crisostomi, R. Gregory and N. Stergioulas, _Rotating Black Holes in Higher Order Gravity_ , _Phys. Rev. D_ 100 (2019) 084020 [1903.05519].
  * (60) K. Takahashi and H. Motohashi, _General Relativity solutions with stealth scalar hair in quadratic higher-order scalar-tensor theories_ , _JCAP_ 06 (2020) 034 [2004.03883].
  * (61) R. C. Bernardo and I. Vega, _Stealth black hole perturbations in kinetic gravity braiding_ , _J. Math. Phys._ 62 (2021) 072501 [2007.06006].
  * (62) M. A. Gorji, H. Motohashi and S. Mukohyama, _Stealth dark energy in scordatura DHOST theory_ , _JCAP_ 03 (2021) 081 [2009.11606].
  * (63) E. Babichev, C. Charmousis, G. Esposito-Farèse and A. Lehébel, _Hamiltonian unboundedness vs stability with an application to Horndeski theory_ , _Phys. Rev. D_ 98 (2018) 104050 [1803.11444].
  * (64) K. Takahashi, H. Motohashi and M. Minamitsuji, _Linear stability analysis of hairy black holes in quadratic degenerate higher-order scalar-tensor theories: Odd-parity perturbations_ , _Phys. Rev. D_ 100 (2019) 024041 [1904.03554].
  * (65) C. de Rham and J. Zhang, _Perturbations of stealth black holes in degenerate higher-order scalar-tensor theories_ , _Phys. Rev. D_ 100 (2019) 124023 [1907.00699].
  * (66) H. Motohashi and S. Mukohyama, _Weakly-coupled stealth solution in scordatura degenerate theory_ , _JCAP_ 01 (2020) 030 [1912.00378].
  * (67) J. Khoury, M. Trodden and S. S. C. Wong, _Existence and instability of hairy black holes in shift-symmetric Horndeski theories_ , _JCAP_ 11 (2020) 044 [2007.01320].
  * (68) K. Tomikawa and T. Kobayashi, _Perturbations and quasinormal modes of black holes with time-dependent scalar hair in shift-symmetric scalar-tensor theories_ , _Phys. Rev. D_ 103 (2021) 084041 [2101.03790].
  * (69) K. Takahashi and H. Motohashi, _Black hole perturbations in DHOST theories: master variables, gradient instability, and strong coupling_ , _JCAP_ 08 (2021) 013 [2106.07128].
  * (70) S. Mukohyama, K. Takahashi and V. Yingcharoenrat, _Generalized Regge-Wheeler equation from Effective Field Theory of black hole perturbations with a timelike scalar profile_ , _JCAP_ 10 (2022) 050 [2208.02943].
  * (71) J. Khoury, T. Noumi, M. Trodden and S. S. C. Wong, _Stability of hairy black holes in shift-symmetric scalar-tensor theories via the effective field theory approach_ , _JCAP_ 04 (2023) 035 [2208.02823].
  * (72) A. De Felice, S. Mukohyama and K. Takahashi, _Approximately stealth black hole in higher-order scalar-tensor theories_ , _JCAP_ 03 (2023) 050 [2212.13031].
  * (73) K. Nakashi, M. Kimura, H. Motohashi and K. Takahashi, _Black hole perturbations in higher-order scalar–tensor theories: initial value problem and dynamical stability_ , _Class. Quant. Grav._ 39 (2022) 175003 [2204.05054].
  * (74) J. Ben Achour, D. Langlois and K. Noui, _Degenerate higher order scalar-tensor theories beyond Horndeski and disformal transformations_ , _Phys. Rev. D_ 93 (2016) 124005 [1602.08398].
  * (75) H. Motohashi, T. Suyama and K. Takahashi, _Fundamental theorem on gauge fixing at the action level_ , _Phys. Rev. D_ 94 (2016) 124021 [1608.00071].
  * (76) S. Mukohyama, K. Takahashi, K. Tomikawa and V. Yingcharoenrat, _Quasinormal modes from EFT of black hole perturbations with timelike scalar profile_ , _JCAP_ 07 (2023) 050 [2304.14304].
  * (77) S. Mukohyama and V. Yingcharoenrat, _Effective field theory of black hole perturbations with timelike scalar profile: formulation_ , _JCAP_ 09 (2022) 010 [2204.00228].
  * (78) R. H. Price, _Nonspherical perturbations of relativistic gravitational collapse. 1. Scalar and gravitational perturbations_ , _Phys. Rev. D_ 5 (1972) 2419.
  * (79) C. Gundlach, R. H. Price and J. Pullin, _Late time behavior of stellar collapse and explosions: 1. Linearized perturbations_ , _Phys. Rev. D_ 49 (1994) 883 [gr-qc/9307009].
  * (80) M. Walker, _Block Diagrams and the Extension of Timelike Two‐Surfaces_ , _Journal of Mathematical Physics_ 11 (2003) 2280.
  * (81) Y. Hatsuda and M. Kimura, _Perturbative quasinormal mode frequencies_ , 2307.16626.