On quasinormal frequencies of blackhole perturbations with external source

We investigate the properties of the quasinormal oscillations in blackhole perturbations with an external source. Though it is not intuitive at a first glimpse, we argue that the source term mainly carries out a similar role as the initial condition. In other words, they do not affect the existing quasinormal frequencies: the characteristic sound of the black holes. It is shown that for various types of perturbations in stationary metrics, the quasinormal frequencies remain unchanged. In particular, for the rotating black hole metric, the derivation involves some subtlety while the same conclusion can still be drawn. Our proof is mostly based on the fact that, under reasonable assumptions, the pole structure of the wave function in the Laplace s-domain remains intact. We claim that it stays valid even when the eigenvalues of the radial and angular parts of the master equation are coupled, as in the case of the Kerr metric. Furthermore, the results are verified numerically by solving the inhomogeneous differential equation and then extracting the dominant complex frequencies by employing the Prony method. We also discuss the scenario where the source term itself is governed by quasinormal oscillations. The implications of the present study to the perturbations in modified gravities, such as the scalar-tensor theories, are addressed.


I. INTRODUCTION
The physical content regarding perturbations in a black hole metric is often viewed analogically as that of a damped harmonic oscillator. Due to the dissipative nature of the system, the frequencies of the oscillations are usually complex. Moreover, when a sinusoidal driving force is applied, the frequency of the steady-state solution is governed by that of the external force, rather than the natural frequency of the oscillator. On the other hand, the problem of black hole quasinormal mode [1][2][3][4][5] is more sophisticated. As an open system, the dissipation demonstrates itself by ingoing waves at the horizon or the outgoing waves at infinity, which subsequently leads to energy loss. Subsequently, for a non-Hermitian system, the relevant excited states are those of quasinormal modes with complex frequencies. Moreover, the boundary condition demands more strenuous efforts, as the solution diverges at both spatial boundaries.
For black hole quasinormal modes, most studies concern the master equation without any external sources. This is largely related to the matter being minimally coupled to the curvature in Einstein's general relativity. However, motivated mainly by the observed accelerated cosmic expansion, theories of modified gravity have become a topic of increasing interest in the last decades. Among many promising possibilities, the latter includes scalartensor [6,7], vector-tensor [8,9], and scalar-vector-tensor [10] theories. There, the matter field can be non-minimally coupled to the curvature sector, and therefore, one might expect the resultant master equations to become inhomogeneous. As an example, the degenerate higher-order scalar-tensor (DHOST) theory [7,11] may admit "stealth" solutions [12]. The latter does not influence the background metric due to its vanishing energy-momentum tensor [13]. The resultant metric differs from the Kerr one by dressing up with a linearly-time dependent field. Indeed, as it has been demonstrated [14,15], under moderate hypotheses, the only non-trivial modification that can be obtained is at the perturbation level. In this regard, the metric perturbations in the DHOST theories have been investigated [16,17] recently. It was shown that the equation of motion for the tensorial perturbations are characterized by some intriguing features. To be specific, the scalar perturbation is shown to be decoupled from those of the Einstein tensor. This result leads to immediate simplifications, namely, the master equations for the tensor perturbations possess the form of linearized Einstein equations supplemented with a source term. The latter, in turn, is governed by the scalar perturbation. Therefore, it is naturally expected that the study of the related quasinormal modes may provide essential information on the stealth scalar hair, as well as the properties of the spacetime of the gravity theory in question.
The present work is motivated by the above intriguing scenarios and aims to study the properties of quasinormal modes with external sources. To be specific, we investigate the case where the right-hand side of the master equation is supplemented with a source term. In comparison to the case of a driven harmonic oscillator, we show that the existing quasinormal frequencies might not be affected when the external source's frequency spectrum is wellbehaved. The results are shown to be valid for not only static Schwarzschild black holes but also rotating ones. For the latter, the eigenvalues of the radial and angular parts of the master equation are coupled. Therefore the arguments based on contour integral have to be revised accordingly.
The rest of the paper is organized as follows. In the following section, as an analogous toy model, we first discuss the solution of a specific driven harmonic oscillator. More rigorous discussions regarding black hole quasinormal models are given in section III in terms of analysis of the pole structure of the relevant Green function in the Laplace s-domain. It is shown that the latter closely follows the same physical intuition obtained by the toy model. The arguments are then further extended to the case of rotating black holes. In section IV, we show that the above results are in agreement with the numerical calculations. We first solve the inhomogeneous differential equation and then employ the Prony method to numerically extract quasinormal frequencies. The latter is compared against the existing results obtained by the conventional method. The last section is dedicated to further discussions and concluding remarks. In particular, we consider a more sophisticated scenario where the source term itself is characterized by quasinormal oscillations. Possible implications regarding applications to theories of modified gravity are addressed.

II. THE QUASINORMAL FREQUENCIES OF A VIBRATING STRING WITH DISSIPATION
As characterized by complex natural frequencies, a damped oscillator is usually employed to illustrate the physical content for the quasinormal oscillations in a dissipation system. Moreover, regarding a dissipative wave equation, a vibrating string subjected to a driven force is an appropriate analogy. To illustrate the main idea, the following derivations concerns a toy model investigated recently by the authors of Ref. [7]. In this model, the wave propagating along the string is governed by the following dissipative wave equation with a source term, namely, where the the string is held fixed at both ends, and therefore the wave function Ψ(t, x), as well as the source S(t, x), satisfies the boundary conditions respectively, where L denotes the length of the string. Here, the relaxation time τ = const. carries out the role of a simplified dissipation mechanism. If the source term vanishes, Ψ is governed by the superposition of the quasinormal oscillations with complex frequencies ω n given by With the presence of the source term, it is straightforward to verify that the formal solution of the wave equation Eq. (1) reads where we have conveniently expanded the source term S(t, x) in the form We note that the integral range for the variable ω is from −∞ to ∞, as it originates from a Fourier transform. Subsequently, for t > 0, according to Jordan's lemma, one chooses the contour to close the integral around the lower half-plane. As a result, the integral evaluates to the summation of the residue of the integrant, namely, Here we have assumed that the source term is a moderately benign function in the sense that S (ω, n) does not contain any singularity. The quasinormal frequencies can be read off by examing the temporal dependence of the above results, namely, the e iω ± n t factors. They are, therefore, determined by the residues of the integrant. The latter is related to the zeros of the denominator on the first line of Eq. (7), which are precisely the frequencies given in Eq. (4).
It is observed that both quasinormal frequencies given by Eq. (4) are below the real axis, and therefore will be taken into account by the residue theorem. If n 2 π 2 τ 2 /L 2 > 1, the two frequencies lie on a horizontal line, symmetric about the imaginary axis. If, on the other hand, n 2 π 2 τ 2 /L 2 < 1, both frequencies lie on the imaginary axis below the origin.
At a first glimpse, the above results seem to contradict our understanding of the driven harmonic oscillators. To be specific, when a sinusoidal external force is applied, the frequency of the steady-state oscillation is known to be identical to that of the driving force. Moreover, the resonance takes place when the magnitude of the driving frequency matches that of the natural frequency of the oscillator. It is noted that the above picture is actually consistent with Eq. (7). If a given frequency characterizes the driven force, it corresponds to the case where S (ω, n) is dominated by a single frequency, namely, S (ω, n) ∼ δ(ω − ω R ) 1 . Subsequently, this will affect the evaluation of residue in Eq. (7). In particular, it corresponds to a pole in the complex plane at ω = ω R − iǫ, where the additional infinitesimal imaginary part iǫ corresponds to a resonance state with infinite half-life. As a result, the steady-state oscillations obtained by the contour integral will be overwhelmed by the contribution from this pole. In other words, a normal mode will govern the late-time behavior of the system, as expected. We note, however, that it is usually not the physical scenario of our interest, where the external force is not alimented by a steady energy source. On the other hand, since we have assumed a non-singular source S (ω, n) with a moderate frequency spectrum, the resonance states of the system are dictated by the quasinormal frequencies given in Eq. (4).
One might argue that the above analogy can only be treated as a toy model when compared with the problem of black hole quasinormal modes. First, the boundary condition of the problem is different: the solution is divergent at both spatial boundaries. Besides, the system is dissipative not due to localized friction but owing to the energy loss from its boundary, namely, the ingoing waves at the horizon and/or the outgoing waves at infinity. As a result, the oscillation frequencies are complex, which can be further traced to the fact that the system is non-Hermitian. In terms of the wave equation, the term concerning the relaxation time τ is replaced by an effective potential. Nonetheless, in the following section, we show that a similar conclusion can also be reached for black hole quasinormal modes. We first discuss static black hole metric, then extend the results to the case of perturbations in rotating black holes.

A. Schwarzschild black hole metric
For a static black hole metric, the perturbation equation of various types of perturbations can be simplified by using the method of separation of variables χ = Ψ(t, r)S(θ)e imϕ . The radial part of the master equation is a second order differential equation [3,4].
where the effective potential V is determined by the given spacetime metric, spins and angular momentum ℓ of the perturbation. For instance, in four-dimensional Schwarzschild or SAdS metric, for massless scalar, electromagnetic and vector gravitational perturbations, it reads where M is the mass of the black hole, and L represents the curvature radius of AdS spacetime so that the Schwarzschild geometry corresponds to L → ∞. The master equation is often conveniently expressed in terms of the tortoise coordinate, x ≡ r * (r) = dr/f . By expanding the external source in terms of spherical harmonics at a given radius, r or x, the resultant radial equation is given by where S(t, x) corresponds to the expansion coefficient for a given harmonics (ℓ, m).
In what follows, one employs the procedure by carrying out the Laplace transform in the time domain [1,2,18] Subsequently, the resultant radial equation in s-domain readŝ where a prime ′ indicates the derivative regarding x, and the source terms on the r.h.s. of the equation consist of S (s, x) and I(s, x). The latter is governed by the initial condition We note that the lower limit of the integrations in Eqs. (12) is "0". Subsequently,f (s, x) and S (s, x) are not able to capture any detail of Ψ(t, x) for t < 0, unless the latter indeed vanish identically in practice. It is apparent that the above equation falls back to that of the sourceless case by taking S(s, x) = 0 [18]. The solution of the inhomogeneous differential equation Eq. (13) can be formally obtained by employing the Green function method. To be specific, where the Green function satisfies It is straightforward to show that where in asymptotically flat spacetimes, which are bounded with ℜs > 0 The wave function thus can be obtained by evaluating the integral where the integral is carried out on a vertical line in the complex plane, where s = ǫ + is 1 with ǫ > 0. Reminiscent of the toy model presented in the previous section, the discrete quasinormal frequencies are again established by evaluating Eq. (19) using the residue theorem. In this case, one employs extended Jordan's lemma to close the contour with a large semicircle to the left of the original integration path [19]. The integration gives rise to the well-known result where s q indicates the poles inside the counter, "other contributions" are referring to those [20][21][22][23] from branch cut on the real axis, essential pole at the origin, and large semicircle. Therefore, putting all pieces together, namely, Eqs. (19), (20), (15), and (17) lead to where the residues are substituted after the last equality. The above results can be rewritten as with where one considers the case where the initial perturbations has a compact support, in other words, it locates in a finite range x 1 < x ′ < x I and the observer is further to the right of it x > x I . The quasinormal frequencies can be extracted from the temporal dependence of the solution, namely, Eq. (23). Since e isqt is the only time-dependent factor, it is dictated by the values of the residues s q . The locations of the poles s q are entirely governed by the Green function Eq. (17), which, in turn, is determined by the zeros of the Wronskian. Therefore, according to the formal solution Eq. (21) or (23), they are irrelevant to c q , where the source S (s, x) is plugged in. As ℜs q < 0 the wave functions diverge at the spatial boundaries, which can be readily seen by substituting s = s q into Eqs. (18), consistent with the results from the Fourier analysis [2,4,24], as mentioned above. It is observed that from Eq. (23), for given initial conditions I(s q , x), one may manipulate the external driving force S (s q , x) so that only one single mode s q presented in the solution.
We note that the above discussions follow closely to those in the literature (see, for instance, Ref. [2,18]). The only difference is that one subtracts the contribution of the external source, namely, S (s, x), from the initial condition I(s, x ′ ) in Eq. (15). It is well-known that the initial conditions of perturbation are irrelevant to the quasinormal frequencies, which characterize the sound of the black hole. In this context, it is inviting to conclude that the external source term on the r.h.s. of the master equation Eq. (11) bears a similar physical content. The Laplace formalism employed in this section facilitates the discussion. On the other hand, from Eq. (23), one finds that the quasinormal modes' amplitudes will still be affected by the external source. Overall, regarding the detection of quasinormal oscillations, the inclusion of external source does imply a significant modification of observables, such as the signal-noise ratio (SNR).
Before closing this subsection, we briefly comment on the equivalence between the above formalism based on Laplace transform and those in terms of Fourier analysis. The results concerning the contour of integral and the quasinormal modes can be compared readily by taking s = −iω [23]. To be more explicit, if one employs the Fourier transform together with the Green function method to solve Eq. (8), the formal solution has the form [25] where one considers the case without a source, and the contributions from the boundary at spatial infinity are irrelevant physically and have been ignored. The Green function is the defined by If we assume that the perturbations vanish identically for t < 0, in other words, G(t, x, x ′ ) = 0 for t < 0. By employing the Fourier transform in the place of Laplace transform, we havẽ whereG Now, it is apparent that, up to an overall sign, the solution of Eq. (27) is essentially identical to Eq. (17) in terms of Laplace transform. The boundary contributions to the formal solution Eq. (26) are precisely those that the initial condition I contribute to Eq. (15). As discussed above, the main reason to employ the Laplace transform is that the formalism provides a transparent interpretation of the role taken by the external source.

B. Kerr black hole metric
As most black holes are likely to be rotating, calculations regarding stationary but rotating metrics are of potential significance from an experimental viewpoint. In this subsection, we extend the above arguments to the case of the Kerr metric. Here, the essential point is that the master equation of the Kerr metric cannot be rewritten in the form of a single second-order ordinary differential equation, such as Eq. (11). To be specific, by employing the method of separation of variables χ = e −iωt e imϕ R(r)S(θ), in standard Boyer-Lindquist coordinates, the master equation is found to be [26] ∆ −s d dr ∆s +1 d dr R (ω, r) + VR(ω, r) = 0, where V (r) = 1 ∆(r) (r 2 + a 2 ) 2 ω 2 − 4Mamωr + a 2 m 2 + 2ia(r − M)ms − 2iM(r 2 − a 2 )sω + (−a 2 ω 2 + 2iωsr −sA ℓm ), ∆(r) = r 2 − 2Mr + a 2 , u ≡ cos θ.
Also, M and aM ≡ J are the mass and angular momentum of the black hole, m ands are the mass and spin of the perturbation field. The solution of the angular part,sS ℓm = s S ℓm (aω, θ, φ), is known as the spin-weighted spheroidal harmonics. Here we have adopted the formalism in Fourier transform for simplicity. Although both equations are ordinary differential equations, the radial equation for the quasinormal frequency ω now depends explicitly onsA ℓm . The latter is determined by the angular part of the master equation, which again involves ω. Therefore, when an external source is introduced, it seems one can no longer straightforwardly employ the arguments presented in the last section. In particular, the arguments based on contour integral seem to work merely for the case where the radial equation is defined in such a way that it is independent of ω. In what follows, however, we elaborate to show that the existing spectrum of quasinormal frequencies remains unchanged. We divide the proof into two parts.
The starting point is to assume that the solution of the homogeneous Eqs. (28)-(29) is already established. First, let us focus on one particular quasinormal frequency ω = ω n,ℓ,m . For a given ω n,ℓ,m , the angular part Eq. (29) is well-defined, and its solution is the spinweighted spheroidal harmonics,sS ℓm . The latter is uniquely associated with a given value ofsA ℓm . Now let us introduce an external source to the perturbation equations (28)- (29). One can show that ω = ω n,ℓ,m must also be a pole of the Green function of the radial part of the resultant master equation. The proof proceeds as follows. It is known that the spinweighted spheroid harmonics form a complete, orthogonal set for a given combination of s, aω, and m [27]. Therefore, it can be employed to expand any arbitrary external source. The expansion coefficient S is a function of radial coordinate r will enter the radial part of the master equation, namely, while the angular part Eq. (29) remains the same. It is note thatsA ℓm is given with respect to given ω n,ℓ,m , thus is denoted bysA ℓmωn . Now, one is allowed to release and vary ω in Eq. (31) in order to solve an equation similar to Eq. (13). As discussed in the last section, one may utilize the Green function method, namely, for real values of ω one solves and then considers analytic continuation of ω onto the complex plane. It is evident that ω n,ℓ,m must be a pole of the above Green function. This is because Eq. (32) does not involve the external source S(ω, r), and therefore the poles must be identical to the quasinormal frequencies of related sourceless scenario. As we have already assumed, the latter, Eqs. (28)- (29), have already be solved and ω n,ℓ,m is one of the quasinormal frequencies. Besides, we note that the other poles of Eq. (32) are irrelevant, since they obviously do not satisfy Eq. (29). Moreover, the forms of Eq. (31) as well as the Green function both change once a different value for ω n,ℓ,m is considered. Secondly, let us consider a given ω but of arbitrary value. Again, the angular part of the master equation Eq. (29) is a well-defined as an eigenvalue problem. Subsequently, its solution, the spin-weighted spheroid harmonics, as a complete, orthogonal set for given s, aω, and m, can be utilized to expand the external source. One finds the following radial equation It is noted that the only difference is thatsA ℓm explicitly depends on ω and it is therefore denoted assA ℓmω . AlthoughsA ℓmω is a function of ω, the above equation is still a second order ordinary differential equation in r. In other words,sA ℓmω can be simply viewed as a constant as long as one is solving the differential equation regarding r. Once more, we will employ the Green function method, where the Green function in question satisfies Now, one is left to observe that the pole at ω = ω n,ℓ,m of the Green function Eq. (32) is also a pole for the Green function Eq. (34). The reason is that the pole ω = ω n,ℓ,m of the Green function Eq. (32) corresponds to one of the zeros of the related Wronskian. The latter is an algebraic (nonlinear) equation for ω. Likewise, the poles of the Green function Eq. (34) also correspond to the zeros of a second Wronskian. The latter is also an algebraic equation except that the constantsA ℓmωn is replaced bysA ℓmω , a function of ω. However, sincē s A ℓmω | ω=ω n,ℓ,m =sA ℓmωn , ω n,ℓ,m must also be a zero of the second Wronskian. We, therefore, complete our proof that quasinormal frequencies ω = ω n,ℓ,m are also the poles for the general problem with the external source.
The above results will be verified in the following section against explicit numerical calculations. Moreover, we note that additional poles, besides those originated from the zeros of the Wronskian, might also be introduced owing to the presence of an external source. One interesting example is that they may come from the "quasi-singularity" of the external source. We will relegate these discussions to the last section of the paper.

IV. NUMERICAL RESULTS
In this section, we demonstrate that the results obtained analytically in previous sections agree with numerical calculations. To be specific, we first solve the inhomogeneous differential equations numerically. Subsequently, the evolution of perturbations in the time domain is used to extract the dominant complex frequencies by utilizing the Prony method. These frequencies are then compared against the numerical results of the corresponding quasinormal modes, obtained by standard approaches.
We first demonstrate the precision of our numerical scheme by studying the toy model presented in section II. Then we proceed to show the results for the Schwarzschild as well as Kerr metrics.  For the toy model, one considers the master Eq. (1) numerically for τ = 1, L = 1, n = 1 and the source Eq. (6) where and Our first goal is to find the temporal dependence of the solution for the two arbitrary sources chosen above. This is accomplished by first solve Eq. (1) in the frequency space and then carry out an inverse Fourier transform at an arbitrary given position x for various time instant t. Although part of the above procedure can be obtained analytically, we have chosen to adopt the numerical approach, since later on, for more complicated scenarios, we will eventually resort to the "brutal" numerical force. The resultant time series are shown in Fig. 1. It is observed that the temporal evolution indeed follows the pattern of quasinormal oscillations.
In order to extract the quasinormal frequencies, the Prony method [28] is employed. The method is a powerful tool in data analysis and signal processing. It can be used to extract the complex frequencies from a regularly spaced time series. The method is implemented by turning a non-linear minimization problem into that of linear least squares in matrix form. As shown below, in practice, even a small dataset of 40 points is often sufficient to extract precise results. In the following, we choose the modified least-squares Prony [28] over others, as the impact of noise is not significant in our study.
For Eq. (35), the two resultant dominating quasinormal frequencies are found to be ω ± (1) = −0.999i − 2.982, −0.999i + 2.967. For Eq. (36), one obtains ω ± (2) = −0.999999i − 2.978190, −0.999998i + 2.978189. When compared with the analytic values ω ± = −i ± √ π 2 − 1 ∼ −i ± 2.978188, one finds that desired precision has been achieved. Next, one proceeds to the case of the Schwarzschild black hole. Here, we consider massless scalar perturbation with the following source term where we takes = 0, r h = 2M = 1, ℓ = 1, L = ∞, V and f are given by Eq. (9)-(10), the tortoise coordinate x = dr/f . It is noted that the factor e iωr V (r)/f 2 (r) is introduced to guarantee that the source satisfies appropriate boundary conditions. The remaining factor 1 1+ω 2 1 r can largely be chosen arbitrarily. To find the temporal evolution, we again solve the master equation in the frequency domain of Eq. (11) by employing a adapted matrix method [29,30]. To be specific, the radius coordinate is transform into a finite interval x ∈ [0, 1] by r → 2M 1−x , which subsequently discretized into 22 spatial grids. For simplicity, we consider α = 1, ℓ = 1. By expressing the function and its derivatives in terms of the function values on the grids, the differential equation is transformed into a system of linear equations represented by a matrix equation. The solution of the equation is then obtained by reverting the matrix, as shown in the left plot of Fig. 2. Subsequently, the inverse Fourier transform is carried out numerically at a given spatial grid x = 5 21 , presented in the right plot of Fig. 2. As an approximation, the numerical integration is only carried for the range ω ∈ [−20, 20], where a necessary precision check has been performed. By employing the Prony method, one can readily extract the most dominate quasinormal frequencies. The resultant value is ω (3) = −0.5847 − 0.1954i, consistent with ω n=0,ℓ=1 = −0.5858 − 0.1953i obtained by the matrix method [30]. Now, we are ready to explore the master equation Eq. (31) for Kerr metric with the following form for the source term where r + = M + √ M 2 − a 2 is the radius of the event horizon. Here, the form r(r−r + ) ∆ e iωr is to guarantee that the external source vanishes at the spatial boundary as a → 0, so that the asymptotical behavior of the wave function remains unchanged. Also, the factor 1 1+ω 2 is again introduced, based on the observation that its presence in Eq. (36) has led to better numerical precision. The latter is probably due to that the resultant numerical integration regarding the inverse Fourier transform converges faster. This choice turns out to be particularly helpful in the present scenario where the numerical precision becomes an impeding factor. In the following calculations, we choose M = 0.5, a = 0.3, and ℓ = 2.
Based on the matrix method, the entire range of the spatial and polar coordinates r and θ is divided by 22 grids. Subsequently, the radial, as well as angular parts of the master equation, are discretized into two matrix equations [31]. We first solve the angular part of the master equation Eq. (29) for a given ω to obtainsA ℓmω . This can be achieved with relatively high precision, namely, with a WorkingPrecision of 100 in Mathematica. The obtained ω to obtainsA ℓmω is substituted back into Eq. (31) to solve for the wave function in the frequency domain. To improve efficiency, we only carry out the calculation for a given spatial point at x = 4 21 , without losing generality. The resultant wave function is shown in the left plot of Fig. 3. To proceed, we evaluate the wave function at 600 discrete points between −30 < ω < 30 and then use those values to approximate the numerical integration in the frequency domain. The resultant time series with 40 points are shown in the right plot of Fig. 3. By using the Prony method, the most dominant quasinormal frequency is found to be ω (4) = −0.9981 − 0.1831i, in good agreement with the value ω n=0,ℓ=2 = 0.9918 − 0.1869i obtained by the 21th order matrix method [31].

V. FURTHER DISCUSSIONS AND CONCLUDING REMARKS
To summarize, in this work, we study the role of possible external sources in blackhole perturbations. We show that even with the presence of the source term, the quasinormal frequencies largely remain unchanged. The statement is valid for various types of perturbation in both static and/or stationary metrics. Although, for rotating black holes, the arguments are elaborated with additional subtlety. The findings are then attested against the numerical calculations for several particular scenarios.
In our discussions, we have ignored a particularly intriguing scenario where the source term itself is governed by quasinormal oscillations. It is understood that a quasinormal frequency with a positive real part implies instability as the oscillations increase exponentially in time. Therefore, one only considers the case that the real part of the s = −iω is negative. Subsequently, it indicates that, in the Laplace s-domain, the source term possesses poles located on the left of the imaginary axis. Since these poles also located inside the contour in Eq. (20), according to the residue theorem, they will introduce additional quasinormal frequencies to the temporal oscillations. The above contributions have not been taken into account in the present work. Also, in our discussions, the effects of the branch cut on the negative real axis have not been considered. However, they are not relevant to the quasinormal frequencies in the context of the present study. The discontinuity from the branch cut arises from that of the solution of the homogeneous radial equation, which satisfies the boundary condition at infinity. It remains unchanged as the external source is introduced. Moreover, as the branch cut stretches from the origin, it primarily associated with the late-time behavior of the perturbations.
Possible implications of the present study include the perturbations in modified gravity theories, such as the scalar-tensor theories. As mentioned above, one relevant feature of the theory is that the scalar perturbations are entirely decoupled from those of the Einstein tensor. Moreover, the metric perturbations in the DHOST theory are found to possess a source term [16,17]. In particular, the master equation for scalar perturbations is shown to be a first-order differential equation decoupled from the Einstein tensor perturbations. As a result, one may obtain the general solution (see, for example, Eq. (26) of Ref. [17]), which does not contain any pole in the frequency domain. In other words, the discussions in section III.B can be readily applied to this case. In this regard, we have demonstrated that while the magnitude of the perturbation wave function will be tailored by the source, the quasinormal frequencies might stay the same. Therefore, the findings of the present work seem to indicate a subtlety in extracting information on the stealth scalar hair in the DHOST theory via quasinormal modes. Further studies along this direction are in progress.