Quasi-Normal Modes of Covariant Loop Quantum Black Holes

Introduction

During the past century, general relativity (GR) has passed a wide range of experimental tests since its birth (see, e.g., [1]). These include the observation of gravitational light deflection during the solar eclipse of 1919, tests of the strong equivalence principle motivated by Nordtvedt’s work in 1968 [2], and measurements of the orbital decay of the binary pulsar PSR B1913+16, discovered by Hulse and Taylor in 1974, to name a few. So far, GR remains one of the most successful theories of gravity because of its predictive accuracy over a wide range of scales. Nevertheless, the nature of dark matter, dark energy, and the quantization of gravity remain open questions. These questions motivate both further tests of GR and the development of theories that go beyond it. The history of modified theories of gravity can be traced back to the early 1920s [2]. Many such theories agree with solar system tests, where gravity is relatively weak, but can differ from GR in other regimes. Observations of binary pulsars have therefore made it possible to rule out or tightly constrain many alternatives [1], while others remain viable and help us explore the nature of gravity. From the Lagrangian, or equivalently the action, point of view, one can modify GR by introducing higher derivatives [3], adding additional fields [4], or relaxing assumptions such as local Lorentz invariance [5, 6].

On the other hand, loop quantum gravity (LQG) provides an important approach to the quantization of gravity. By constructing a quantum theory of Riemannian geometry, it offers a framework for investigating quantum gravitational effects, including the possible resolution of classical singularities [7]. Ideas and methods from LQG have also been applied to symmetry-reduced cosmological models in loop quantum cosmology (LQC) [8, 9, 10]. This has enabled detailed studies of quantum corrections and singularity resolution in these models [8, 9, 10, 11]. In widely studied homogeneous LQC models, quantum geometrical effects replace the classical big bang or big crunch singularity by a regular bounce, connecting contracting and expanding branches of spacetime.

Black holes (BHs) are among the most mysterious objects in the universe and are a major prediction of GR [12]. Interest in BHs has grown further with the Event Horizon Telescope (EHT) observations of the shadows of M87* and Sagittarius $A^{*}$, the latter at the center of our own Milky Way galaxy [13, 14]. Many galaxies are believed to host supermassive BHs at their centers [15, 16]. Their strong gravitational fields and interactions with surrounding matter provide opportunities to study phenomena such as relativistic jets and tidal disruption. From the theoretical point of view, BHs and other strong-field systems also provide opportunities for testing theories of gravity beyond GR (see, e.g., [17, 18, 19, 20, 21, 22, 23, 24]). This role also extends to investigating possible quantum gravitational effects in effective BH models motivated by LQG [25, 26]. Under suitable assumptions, the classical singularity theorems imply geodesic incompleteness, indicating limitations of the classical description [27]. Studies of effective LQG models suggest that quantum geometry can resolve the spacelike singularity inside BHs. Singularity resolution has also been extensively studied in LQC [8, 9, 10]. Such studies often use effective equations to describe how quantum corrections modify the classical geometry.

The test of gravitational theories continues in the strong-field regime. In particular, the first direct detection of gravitational waves (GWs) from the coalescence of two BHs by Advanced LIGO marked the beginning of GW astronomy [28]. Since then, the LIGO/Virgo/KAGRA (LVK) collaborations have reported hundreds of compact binary coalescence candidates, with catalogs and public data supporting increasingly detailed studies [29, 30, 31, 32, 33, 34, 35]. Planned ground-based detectors and space-based missions aim to improve sensitivity and extend the frequency ranges accessible to GW observations [36, 37]. These include the Cosmic Explorer [38], Einstein Telescope [39], LISA [40], TianQin [41, 42, 43], Taiji [44], and DECIGO [45]. Together, they could probe a wider range of sources and distances. These prospects have stimulated studies of BH quasinormal modes (QNMs) [46, 47], extreme mass ratio inspirals (EMRIs) [48, 49, 50], and the effects of dark matter (DM) on BH dynamics and GW signals [16, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], to name a few.

During the ringdown stage of a coalescence, the newly formed BH approaches a stationary state. When linear perturbation theory applies, its response can be approximately described by a superposition of QNMs [61, 62, 63, 64, 65, 66], although nonlinear contributions can also be important, especially in the early ringdown [67, 68]. Scalar-field, vector (electromagnetic), and tensor (gravitational) perturbations and their associated QNM problems have been extensively studied in GR and in modified theories of gravity. In particular, test scalar fields provide a relatively simple setting for studying how the background geometry affects wave propagation and QNM spectra [46]. From the theoretical point of view, QNMs are characteristic modes of dissipative systems with generally complex frequencies $\omega$. The real part determines the angular frequency of oscillation, while the imaginary part determines the decay or growth rate. For a damped mode with time dependence $e^{-i\omega t}$, the damping time is $1/\lvert\operatorname{Im}\omega\rvert$. These complex frequencies carry information about the background geometry and perturbation dynamics, providing a way to test gravity in the strong-field regime. In four-dimensional GR, the standard no-hair result states that stationary, asymptotically flat electrovacuum BHs are characterized by mass, charge, and angular momentum under suitable assumptions [46, 47]. For a specified perturbing field and mode, their QNM frequencies are therefore fixed by these parameters. Measuring multiple modes can test the consistency of these relations. The dependence of QNMs on the background geometry and perturbation sector has also been explored in specific BH models [16, 51, 69]. An important example is the isospectrality of axial and polar gravitational perturbations of the Schwarzschild BH [66]. For more details, we refer readers to [46, 47, 70, 71].

On the other hand, from the observational point of view, identifying modes beyond the dominant ringdown contribution remains challenging with current detectors [72]. Searches for overtones in GW150914 have led to extensive discussion of the effects of noise, the choice of fitting interval, and the waveform model [73, 74, 75, 76, 77, 78]. Evidence for a subdominant angular mode has also been reported for GW190521 [79]. The high signal-to-noise ratio of GW250114 has provided stronger evidence for the first overtone of the dominant quadrupolar mode and enabled more precise tests of the Kerr spectrum [80, 81, 82]. Future third-generation ground-based detectors and space-based missions are expected to improve the prospects for resolving subleading modes in sufficiently loud events [47, 83]. Interpreting such observations within theories beyond GR requires predictions for how changes in the background geometry and perturbation dynamics affect BH ringing. We study this question in an effective LQG framework using scalar perturbations as a first step. Although these perturbations do not directly describe GWs, they provide a simpler setting in which to examine how holonomy corrections affect the QNM spectrum and the time-domain response. Such studies can probe the dynamical consequences of the effective theory’s holonomy modifications [84]. This study forms part of a broader research program that we aim to extend to gravitational perturbations and their observable GW signatures.

This paper focuses on a class of covariant loop quantum black hole solutions [85, 86], obtained by applying emergent modified gravity to holonomy modifications motivated by LQG [87, 88]. To distinguish them from other loop quantum black holes (LQBHs), we refer to them as cLQBHs. In the class considered here, the solutions are characterized by two functions $\chi(r)$ and $\lambda(r)$. We restrict attention to $\chi(r)=1$ and focus on the localized profile of the dimensionless function $\lambda(r)$ in case (ii), with the previously studied decreasing profile in case (i) included for reference,

where $M>0$ is the mass parameter, $\Delta>0$, $\epsilon>0$, $\sigma>0$, and $r_0 > 2M$ (with $c = G = 1$).¹

By considering a massless, minimally coupled test scalar field on these fixed backgrounds [84], we investigate the corresponding wave dynamics and QNM spectra as a first step toward future studies of gravitational perturbations. The decreasing profile in case (i) provides a reference from previous discussions of covariant holonomy modifications [93]. Our focus is on case (ii), whose Gaussian factor allows us to vary the position, width, and strength of the resulting deformation of the effective potential, as we demonstrate below. The same background geometry determines both the effective potential and the relation between the areal radius $r$ and the tortoise coordinate. This enables us to examine how localized holonomy corrections affect BH ringing, with a view to future investigations of their possible observational implications [94, 95].

It is important to note that the localized correction in case (ii) produces a bump structure in the scalar effective potential, connecting this model to the study of the bump effect on BH QNMs [90]. Recent work in this direction has explored spectral sensitivity, mode excitation, and the relation between frequency-domain spectra and time-domain waveforms [91, 96, 97, 98, 99, 100, 101, 102]. These studies have identified several characteristic features. In particular, a small bump far from the main potential peak can strongly shift QNM frequencies while leaving the early ringdown waveform nearly unchanged, and its effects on the waveform may become more apparent at later times [91, 97]. Related studies have also examined whether localized negative or time-dependent potential deformations can trigger dynamical instability [103]. Some of these studies use toy models with prescribed potential bumps or simplified scattering problems. Meanwhile, other studies

examine complex potential structures, including double peaks in specific gravitational theories [104, 105], as well as modifications produced by modeled matter environments [106, 107]. These approaches provide complementary insights into the origin and consequences of localized potential structures. In this context, our chosen case (ii) provides a concrete cLQBH realization of the bump effect. As noted above, the bump follows from the cLQBH geometry specified by the holonomy function $\lambda(r)$, rather than being imposed directly on the GR effective potential. This gives the bump a definite origin within an effective gravitational theory and allows systematic calculations of both the QNM spectrum and the time-domain response. It thus provides a bridge between the bump effect, BH spectroscopy, and observational tests of gravity. The delayed effects discussed above also motivate a quantitative study of when the waveform departs from its GR counterpart, using a suitable deviation criterion. The physical significance of this timescale also depends on whether the deviation occurs while the signal still has appreciable amplitude. Considering the deviation time together with the complex QNM frequencies $\omega$ and the signal decay helps assess the observational significance of the bump effect, including exploratory comparisons of frequency and damping scales with published GW measurements. These considerations motivate the study of an “observation window” in which deviations may become apparent before the signal becomes too weak, bringing the bump effect into a concrete framework for QNM calculations and tests of gravity.

The rest of the paper is organized as follows. Section II briefly introduces the cLQBH solutions and the master equation for a test scalar field, along with the corresponding effective potential. In Sec. III, we study QNMs in both the frequency and time domains [mainly for case (ii), with case (i) as a reference] and compare the results with their Schwarzschild counterparts. We also discuss their dependence on the model parameters and make preliminary comparisons with published LIGO data to explore possible constraints on the model. In Sec. IV, we discuss the deviation time and its relation to waveform decay and consider their observational relevance, including the idea of an observation window and the role of detector properties. Finally, our conclusions and outlook are given in Sec. V.

Before proceeding further, we note that the scalar perturbations studied here do not directly describe the spin-2 GWs emitted during BH ringdown [46, 47, 70, 71]. Our main comparison is therefore between the same test scalar field on cLQBH and Schwarzschild backgrounds. Comparisons with published GW data are exploratory, while quantitative constraints from GW ringdown require gravitational perturbation equations consistent with the cLQBH dynamics and a statistical analysis that accounts for detector noise (see, e.g., [47, 80]). For future extensions to gravitational perturbations, useful technical references include the master equations and QNM calculations for odd-parity perturbations in other effective LQBH models [108, 109]. Work covering both odd- and even-parity perturbations has also revealed a breaking of isospectrality [110]. Gravitational perturbation theory has also been developed at the quantum level through canonical formulations and hybrid quantization of both parity sectors on Kantowski–Sachs backgrounds relevant to BH interiors [111, 112]. More recently, this framework has been extended to BH exteriors [113, 114, 115] and applied to gravitational QNM calculations in an effective LQBH model [116]. More directly related to our scalar analysis, Bojowald et al. [84] derived scalar master equations within the cLQBH framework and compared minimally and nonminimally coupled QNMs for constant and decreasing holonomy functions using a third-order WKB approximation. Their results show that both the holonomy profile and the matter coupling can substantially affect the QNM spectrum. Their minimally coupled equations provide a basis for our analysis, with case (i) retained as the reference discussed above. Building on this work, we focus on the Gaussian-localized holonomy profile in case (ii). We study the resulting bump effect in both the frequency and time domains, linking spectral changes, waveform deviations, and their possible observational significance.

Throughout this paper, we adopt units with $c=G=1$ [84], where $c$ is the speed of light and $G$ is the gravitational constant. Where convenient, we further set $M=1$, completing the choice of units for length, time, and mass. All dimensional quantities must be expressed consistently in this unit system. Greek indices run from 0 to 3. Other index conventions will be specified where needed. We adopt the metric signature $(-,+,+,+)$.

Spherically symmetric cLQBHs and scalar field perturbations

In this section, we first briefly introduce the static, spherically symmetric cLQBH solutions constructed in [85], and then consider linear perturbations of a massless, minimally coupled test scalar field on these backgrounds. After rescaling the time coordinate by a constant factor, the line element can be written as

where $r$ is the areal radius, $d\Omega^2\equiv d\theta^2+\sin^2\theta\,d\varphi^2$, and

Here $\chi(r)$ and $\lambda(r)$ are two real functions of $r$, with $\chi(r)\ne0$ in the region under consideration. The constants $M$ and $\Lambda$ denote the mass parameter and cosmological constant, respectively. In Ref. [86], Hawking radiation was studied for the decreasing function

[cf., case (i) in (1.1)], with $\Lambda=0$ and $\chi(r)=1$, where $\Delta>0$ is a constant with dimensions of area related to the area scale in LQG. For these choices of $\Lambda$ and $\chi(r)$, $g_{tt}=-f(r)$ coincides with the corresponding Schwarzschild metric component [12], and the metric horizon (MH) is located at $r=r_{\mathrm{MH}}\equiv2M$ for both profiles in (1.1). The covariant framework allows considerable freedom in the choice of $\lambda(r)$, subject to regularity and the required asymptotic behavior. More specific quantum-gravity models may impose additional restrictions [84, 85]. Our choice in case (ii) satisfies these geometric requirements in the exterior and produces a bump in the scalar effective potential, as discussed in the Introduction.

Bojowald et al. [84] studied scalar QNMs for both constant and decreasing functions $\lambda(r)$, considering minimal and nonminimal couplings. The decreasing function used in their minimally coupled analysis corresponds to case (i) in (1.1), and their results provide an important benchmark for our work. However, our main focus is case (ii), in which $\lambda(r)$ has a Gaussian-like radial profile. For the remainder of this paper, we restrict attention to

Then, following [94], we consider scalar perturbations governed by the Klein-Gordon equation

where $\Phi=\Phi(t,r,\theta,\varphi)$ denotes the scalar field in Boyer–Lindquist-type coordinates [117] and $g$ is the determinant of the background metric. This equation is compatible with the covariance of the theory [84, 86]. Separating variables as $\Phi=Y_{lm}(\theta,\varphi)\Psi(t,r)/r$, where $Y_{lm}$ are spherical harmonics, we obtain the master equation for $\Psi$,

where the effective potential $V_{\mathrm{eff}}$, with dimensions of $M^{-2}$, is given by [118]

where $L=l(l+1)$, $V_{\mathrm{GR}}=(r-2M)(Lr+2M)/r^4$ denotes the effective potential for a massless scalar field on the Schwarzschild background [84]², and $x$ is the tortoise coordinate. coordinate defined through

Here a prime denotes differentiation with respect to $r$. This effective potential agrees with that given in [86] after accounting for the different notation. Since the background is spherically symmetric, the radial equation and QNM frequencies are independent of the angular index $m$ [66]. Thus, without loss of generality, we set $m=0$. For given functions $f(r)$ and $\lambda(r)$ and a chosen value of $l$, we solve (6) for $\Psi$.

Notice that, the master equation (6) mainly serve the time-domain calculations for the QNMs. For the frequency-domain analysis the calculations of $\omega$, a more straightforward way is to use its frequency-domain counterpart by doing the Fourier transformation, with which we obtain

where $\Psi(t,x)=e^{-i\omega t}\psi(x)$. Here, $\omega$ is the complex QNM frequency we are looking for.

For later convenience, and following the notation used in, e.g., [90, 97], let us introduce the “bump potential”

so that $V_{\mathrm{eff}}=V_{\mathrm{GR}}+V_{\mathrm{bump}}$ (notice that no additional overall amplitude factor is placed in front of $V_{\mathrm{bump}}$). Substituting the expressions for $\lambda(r)$ in case (ii) and $f(r)$, we obtain

where, following [94], we have introduced the dimensionless quantities $\tilde r\equiv r/M$ and $\tilde r_0\equiv r_0/M$ to simplify the last line. In fact, using $M$ as the reference scale, the physical quantities of interest considered here can be expressed in dimensionless form, as illustrated by $r$ and $r_0$. For example, according to (11), we can further introduce the dimensionless bump potential $\tilde V_{\mathrm{bump}}\equiv M^2V_{\mathrm{bump}}$. The numerical values of these dimensionless quantities are equivalent to those obtained by setting $M=1$ in the corresponding dimensional expressions. Taking advantage of this equivalence and the simplicity of dimensionless quantities, we shall continue to use the “tilde” notation to construct dimensionless expressions, or set $M=1$, whenever convenient and without causing confusion. Notice that $V_{\mathrm{bump}}$ is independent of $l$. Also notice that, since $f^3\lambda^2$ tends to zero at both boundaries, $r\to r_{\mathrm{MH}}$ and $r\to+\infty$, while being positive at intermediate radii, it follows from (10) that $V_{\mathrm{bump}}$ must contain both positive and negative parts, as shown in the figure to be given below.

We plot a representative case of $M^2V_{\mathrm{bump}}$ in Figure 1 by setting $r_0/M=12$, $\epsilon=0.5$, and $\sigma=0.9$ ($V_{\mathrm{bump}}$ is independent of $l$, as mentioned earlier). This figure illustrates the typical behavior of the bump potential in case (ii) [cf. 1.1]. In particular, the bump exhibits a “sine-like” behavior around $r=r_0$, although the height of the local peak generally differs from the depth of the local trough. In general, the bump also has a nonzero integral, $\int_{-\infty}^{+\infty}V_{\mathrm{bump}}\,dx\ne0$. More importantly, with the assistance of Figure 1, we can define the location, width and amplitude of such bumps.

More specifically, the first line of (11) shows that $V_{\mathrm{bump}}(r=r_{\mathrm{MH}})=V_{\mathrm{bump}}(r=r_0)=0$. The zero crossing at $r=r_0$ naturally allows us to define this point as the center (location) of the bump, as marked by the green arrow in Figure 1. In contrast, the bump’s width and amplitude are not directly specified by $\sigma$ and $\epsilon$, respectively. Although $\epsilon^2$ independently scales the amplitude, changing $\sigma$ generally affects both the width and amplitude. Therefore, as illustrated in Figure 1 (see the red and orange markers), we introduce the dimensionless quantities $\mathcal{A}$ and $\mathcal{W}$ to describe the amplitude and width directly. Together with a given $r_0/M$, these quantities allow us to solve, usually numerically, for the corresponding $\sigma$ and $\epsilon$ using (11). For the cases considered here, this procedure establishes a correspondence between $\{r_0/M,\epsilon,\sigma\}$ and $\{r_0/M,\mathcal{A},\mathcal{W}\}$ through the exterior local peak and trough. This enables us to investigate the individual effects of the location, amplitude and width of the bump, rather than their combined effects when varying the original parameters.

By examining the structure of $V_{\mathrm{bump}}$ in the last line of (11), we notice that, for a given set of $\{r_0/M,\mathcal{A},\mathcal{W}\}$, we can first determine $\sigma$ from $r_0/M$ and $\mathcal{W}$. The same pair of $\{r_0,\mathcal{W}\}$ also determines $\epsilon^2/\mathcal{A}$, whose reciprocal gives the bump amplitude for $\epsilon=1$. The amplitude information is therefore already contained in this quantity, with $\epsilon^2$ providing the remaining overall adjustment. We plot these two dimensionless quantities as functions of $\mathcal{W}$ for several choices of $r_0/M$ in the top and bottom panels of Figure 2, respectively. Once $\mathcal{A}$ is prescribed, obtaining $\epsilon$ from the bottom panel is a trivial algebraic step, so separate numerical plots of $\epsilon$ would add no information to this parameter mapping and are omitted here. This simplification illustrates a broader consideration in our study. The large parameter space entails substantial computational effort and can obscure the principal physical trends among numerous parameter combinations. We therefore organize the investigation around the mathematical structure and the relative importance of the parameter effects, seeking a manageable set of cases that reveals the essential behavior clearly. This consideration will also guide the subsequent analysis. Returning to Figure 2, we observe that both $\sigma$ and $\epsilon^2/\mathcal{A}$ depend approximately linearly on $\mathcal{W}$ for a given $r_0/M$, while $\sigma$ is insensitive to $r_0/M$.

We now return to the features of $V_{\mathrm{bump}}$ and $V_{\mathrm{eff}}$. Motivated by the observational significance of the $l=2$ gravitational-wave multipole [119], and following the conventions of, e.g., [90, 91, 94], our analysis will mainly focus on the representative $l=2$ case in the following. Figure 3 shows the behavior of the dimensionless potentials $M^2V_{\mathrm{bump}}$ and $M^2V_{\mathrm{eff}}$ for 8 carefully selected typical combinations of $\{r_0,\mathcal{A},\mathcal{W}\}$ at $l=2$ (The factor $M^2$ was multiplied to $V_{\mathrm{eff}}$ to balance its original dimension so that the overall $\widetilde{V}_{\mathrm{eff}}=M^2V_{\mathrm{eff}}$ can be dimensionless). These

choices will serve as starting points for subsequent explorations of the parameter space spanned by $\{r_0,\mathcal{A},\mathcal{W}\}$. In particular, as illustrated in Figure 3, this paper mainly considers cases where $V_{\mathrm{eff}}(r)>0$ throughout the exterior region $r\in(2M,+\infty)$. This is a reasonable choice, since negative regions of the effective potential can favor the development of dynamical instabilities [103, 120, 121], an issue beyond the scope of the present work (although a locally negative effective potential does not necessarily imply dynamical instability [122], and this restriction does not exclude a locally negative $V_{\mathrm{bump}}$ from our study).

Before proceeding to the QNM calculations and the associated analysis, let us briefly examine the behavior of the tortoise coordinate [cf. (8)]. Substituting the expressions for $f(r)$ and $H(r)$, we obtain

where the last line gives the dimensionless form, with $\widetilde{x}\equiv $x/M$ and $dr/dx=d\widetilde{r}/d\widetilde{x}$. For the Schwarzschild background, $dr/dx=1-2M/r$, and a convenient choice of the integration constant gives $x_{\mathrm{GR}}=r+2M\ln[(r-2M)/M]$ [66, 70] so that $r=3M$ leads to $x_{\mathrm{GR}}=3M$.³ We can integrate (12) numerically to obtain the relation between $r$ and $x$ for a given $\{r_{0},\epsilon,\sigma\}$, or equivalently, $\{r_{0},\mathcal{A},\mathcal{W}\}$. Figure 4 shows the typical behavior of $x(r)$ for $\{r_{0}/M,\mathcal{A},\mathcal{W}\}=\{12,0.07,0.5\}$. To compare the cLQBH and GR coordinates, we use the freedom in the integration constant to impose $\widetilde{x}(\widetilde{r}=3)=3$ for both backgrounds, consistently with the expression for $x_{\mathrm{GR}}$ above. This merely fixes the origin of $x$, without affecting the physical results [51]. Clearly, the $x$ maps

$r/M\in(2,+\infty)$ to $x/M\in(-\infty,+\infty)$. In the bottom panel of Figure 4, we further plot the dimensionless difference $\Delta x/M\equiv(x_{\mathrm{cLQBH}}-x_{\mathrm{GR}})/M$ as a function of $r/M$. Together, the two panels show that the bump has a small influence on $x$ in this example, while $\Delta x$ remains approximately constant away from the bump region on either side of $r=r_{0}$. The location $r=r_{0}$ marks the center of the transition between these two nearly constant values. Furthermore, Figure 4 illustrates the mapping between $r$ and $x$, which will be used in the subsequent QNM calculations.

Quasi-Normal Modes of cLQBHs

At this point, we are at the stage of calculating the corresponding QNMs from the master equations given above. Specifically, we can execute the time-domain calculations by starting from (6), and the frequencydomain ones by starting from its counterpart, (9). In principle, for any given $l$ and $\{r_0,A,\mathcal{W}\}$, by integrating (6) we can obtain the temporal evolution of $\Psi(t,x_{\rm obs})$ (where $x_{\rm obs}$ is the location of the observer to be assigned below) with suitable physical boundary conditions to be described at appropriate places. Similarly, for the frequency-domain analysis, one can further extract a set of $\omega$ (the full notation is $\omega_{ln}$, with $l,n=0,1,2,3,\ldots$ and no $m$ index since the background is non-rotating, though we in general omit these subscripts and simply write $\omega$ [46, 84]), either from $\Psi(t,x_{\rm obs})$ or by directly solving (9).

In the literature, various methods are available for handling such QNM problems. Examples include the WKB approximation [63, 65, 123, 124] (with the Padé improvement [125, 126]), the continued fraction method (CFM) [64, 71, 127], the matrix method [51, 128, 129], the pseudospectral method [130, 131], the shooting method [132], and the finite difference method (FDM) [94, 133, 134, 135], together with the Runge–Kutta method [97, 136]. In what follows, we adopt the ones suitable for our calculations in practice. It is worth noting that, while testing these methods, we performed some preliminary calculations and compared them with the case (i) results of [84]. These comparisons served to test our numerical implementation and to improve our understanding of the scalar-field dynamics in cLQBHs. Together with checks against selected results from the studies of the bump effect [97], they guided the development of the numerical codes and the computational workflow employed below.

Frequency-domain analysis

Let us first pay attention to the frequency-domain calculations and relevant analysis for the cLQBH mentioned above. As usual, we shall treat the $l=2$ case as the representative one. However, the selection of the suitable physical parameters $\{r_0,\mathcal{A},\mathcal{W}\}$ can be a more complicated task. First of all, as mentioned earlier, the current work is focusing on the regime where $V_{\rm eff}(r)>0$. This leads to some implicit constraints on $\{r_0,\mathcal{A},\mathcal{W}\}$. By referring to the cases considered in Figure 2 and Figure 3, we shall mainly focus on some subspaces of the parameters that satisfy $\tilde{r}_0\in(5,50)$, $\mathcal{A}\in(0,0.1)$ and $\mathcal{W}\in(0,5)$.

We first notice that we have $V_{\rm eff}(x)\to0$ at the two boundaries of $x$ [cf., (7)], viz., $x\to-\infty$ ($\tilde{r}\to2^+$) and $x\to+\infty$ ($\tilde{r}\to+\infty$). According to (9), that brings us $\psi\propto e^{\pm i\omega x}$ at the two boundaries. The physical boundary conditions, i.e., the pure-ingoing-wave condition at the metric horizon and the pure-outgoing-wave condition at the spatial infinity, can be translated as the following Robin boundary conditions [84, 97, 137]

Therefore, we can pick out the correct branch of $\psi(x)$ as the boundary conditions according to (13) in the following calculations.

The shooting method [132, 138, 139, 140] mentioned above is well suited to our mathematical problem. Choosing this method can make the calculations for $\omega$ straightforward, clean and manageable. Following [132], we introduce $\omega=\omega_R+i\omega_I$ and $\phi=\phi_R+i\phi_I$, where $\psi(x)=e^{i\int\phi(x)\mathrm{d}x}$. Substituting these quantities into (9) leads to the Riccati equations [71]

Here, the subscripts $R$ and $I$ denote the real and imaginary parts, respectively. The coupled complex ordinary differential equations (ODEs) given by (14) are the equations we are going to integrate in our Mathematica code (and its Python counterpart) in practice. To fulfil the integrations, let us construct the numerical boundary values for the near-horizon as well as the far-zone regimes by using (13) and the Frobenius method [141].

In practice, to reduce our computational labor, as mentioned earlier, it is very convenient to work with the dimensionless form of the ODEs or, equivalently, by setting $M=1$. That is, we shall work on top of the following equation

To be specific, we expand $\psi$ near the $\tilde{r}\to+\infty$ regime by

where $\xi\equiv M/r$ and the dimensionless factors $\alpha_k$ will be determined later. Notice that here we have $\omega x=(M\omega)(x/M)=\tilde{\omega}\tilde{x}$. Similarly, we expand $\psi$ near the $\tilde{r}\to2^+$ regime by

where $\zeta\equiv(r-r_{\mathrm{MH}})/M$ and the dimensionless factors $\beta_k$ will be determined later. It is worth mentioning here that one can never expand the polynomials in (16) or (17) to the infinite powers in reality. Instead, we shall expand them to a certain finite order by measuring the tolerance of the numerical error for the calculations of $\omega$. At the same time, due to the relation between $\phi$ and $\psi$, we have the freedom to normalize (16) and (17) to meet our purposes by properly choosing integration constants. Thus, without loss of generality, let us set $\alpha_0=\beta_0=1$.

After that, we insert the polynomials in (16) and (17) into the ODE (15). By solving that order by order, we can arrive at some recursion relations of the coefficients, which read

and

where

with $k\in\mathbb{N}^+$ and $p(r)=p(\tilde{r})\equiv d\tilde{r}/d\tilde{x}=dr/dx$ [cf., (8)]. Notice that the factors $\bar{A}_j$, $\bar{B}_j$ and $\bar{C}_j$ here are nothing but the coefficients for the Taylor expansions of the corresponding functions.

With the above in hand, for any given $l$ and $\{r_0,\mathcal{A},\mathcal{W}\}$, we can generate the numerical value of $\psi(x=x_-)$ [so that $\phi(x_-)$] by using (17) and integrate the ODEs (14) from a chosen $x_-$ [with $r_-=r(x=x_-)$ close to the radius of the metric horizon] to the matching point located at, say, $x=x_m$, at which we obtain $\phi_-\equiv\phi(x=x_m)$. Similarly, we can generate the numerical value of $\psi(x=x_+)$ [so that $\phi(x_+)$] by using (16) and integrate the ODEs (14) from a chosen $x_+$ [with $r_+=r(x=x_+)$ to be sufficiently large] to the same matching point located at $x=x_m$ (a typical choice of $x_m$ is to put its corresponding $r_m$ near the radial position of the main peak of $V_{\mathrm{eff}}$), at which we obtain $\phi_+\equiv\phi(x=x_m)$.⁴

For each iteration, we pick up a reasonable test value of $\omega$ for the numerical integration and check the condition $|\phi_+-\phi_-|\lesssim\epsilon_0$ (which is equivalent to using the Wronskian condition [141, 142]), where $\epsilon_0$ is a small number to be selected to meet the accuracy requirement of our calculations. We iteratively adjust the trial value of $\omega$ until the condition $|\phi_+-\phi_-|\lesssim\epsilon_0$ is met (in practice, a typical choice is $\epsilon_0=2\times10^{-9}$ in units with $M=1$). At that stage, the latest test $\omega$ will be the QNM frequency we are looking for.

At the beginning part of this subsection we have mentioned roughly the parameter space under consideration for this current work (although we are not going to explore every corner of the space mentioned in there). That is to cover two anchor points of the full parameter space, i.e., the one in the vicinity of the GR limit and the one far beyond that (while keeping the deviation as a finite number, as an extremely huge discrepancy is meaningless for our purposes). Here are some considerations regarding the chosen parameter space.

As mentioned earlier, the case (i) [cf., (1)] considered in [84] represents a crucial benchmark for our current study. To do a more direct comparison and by taking the same manner [cf., (7)], we may first construct the counterpart of our $V_{\mathrm{bump}}$ for case (i) by introducing a $\delta V_i\equiv V_i-V_{\mathrm{GR}}$, where $V_i$ denotes the effective potential for the case (i). For the physical parameter $\Delta$, the authors of [84] considered the case where $\Delta/(4M^2)\in[0,1]$. That brings the amplitude of $M^2\delta V_i$ to be able to reach the order of $\mathcal{O}(10^{-3})$. However, they actually did not set up an upper limit for this factor. If, by considering the ultralight BHs [143, 144, 145, 146], the amplitude of $M^2\delta V_i$ can in principle reach the order of approximately $\mathcal{O}(10^{-1})$. This gives us the reason to mainly focus on a sufficiently large $\mathcal{A}$ while keeping it no more than the order of $\mathcal{O}(10^{-1})$, which can bring us to the balanced point between finding an observable deviation from that of GR (at least for some theoretical purposes) and avoiding taking into account unnecessary extreme cases.

On the other hand, since the bump effect on QNM has been investigated intensively during the past several years [90, 91, 97, 101, 147], we shall refer to works for some clues in determining a reasonable range for $r_0$ that is worth exploring. These studies imply that a typical choice of $r_0/M$ can be of the order $\mathcal{O}(10)-\mathcal{O}(100)$. However, we shall not attach equal significance to these choices of $r_0$, since this current work is not a comprehensive study of the bump effect nor purely an extension in this direction. Under the requirement of $V_{\mathrm{eff}}>0$, the allowed range for $\mathcal{A}$ will be too small for a large $r_0$, which will dilute the contribution of the corresponding $\mathcal{A}$ to certain degree. Therefore, for our purpose, as seen from Figure 3, we will mainly focus on the values of $r_0/M$ smaller than 30. In fact, our choices can make some special contributions to the area of the bump effect, as most of the existing works have not paid enough attention to the small $r_0$ regime. Our chosen range for $r_0$ can link the near-zone cases to some of the typical intermediate $r_0$’s existing in the literature. Similarly, we can find the counterpart of our $\mathcal{W}$ in the literature. A typical size for it ranges approximately from 1 to 3. Nonetheless, to extend our study to the “sharp peak” regime, we will allow $\mathcal{W}$ to be smaller than 1 [102, 148].

Knowing all of these, we plot out the real and imaginary parts of the dimensionless quantity $M\omega$ (by setting $l=2$ and $n=0$) for four representative parameter choices, i.e., $\{r_0/M,\mathcal{W}\}=\{6,0.2\},\{16.8,0.2\},\{6,3\},\{16.8,3\}$, as a function of the amplitude $\mathcal{A}$ in Figure 5. In this way, we cover four typical scenarios, namely, a sharp/wide bump located in the near/intermediate zone. By observing Figure 5 for the four independent cases, we notice that the amount of deviation from GR for $\omega$ in general increases monotonically with $\mathcal{A}$, which agrees with our intuition. Similarly, a wider as well as farther bump tends to drag $M\omega$ away from its GR counterpart. Indeed, we also notice that, in the relatively far zone (e.g., $r_0/M\gg2$), the behavior of $M\omega$ can be more complicated. That is partially due to the fact that a relatively large $r_0$ is strongly connected to some fancy behaviors of $\omega$, as revealed by the existing literature [90, 97, 98, 100] and Figure 6. Our results from Figure 6 show that this sort of spiral pattern exists not only for positive bumps (see, e.g., [149]) but also for mixed ones. In contrast, the trajectory of $\omega$ will tend to be more regular with a changing $\mathcal{W}$, as seen from Figure 7. As a connection to some of the previous works (including ours) [138, 150], we would also like to point out another issue arising from Figure 6. Replacing the purely ingoing condition with a partially reflective boundary near the horizon, characterized by a generally complex reflection coefficient, provides another way to modify QNM trajectories [151].

Motivated also by effective descriptions of localized potential perturbations in terms of reflectivity [99], it would be interesting to investigate whether suitable variations of the reflection amplitude and phase can reproduce spiral trajectories resembling Figure 6, and whether a correspondence with the bump parameters can be established, (especially when some sort of degeneracy can occur). We leave this possible connection between potential deformations and modified boundary conditions for future investigation.

It is worth mentioning here that Figure 5 compares our results with approximate reference ranges motivated by LIGO observations and TianQin forecasts. The LVK Collaboration recently reported constraints on departures from the Kerr QNM spectrum using the event GW250114 [80]. Although our calculations concern scalar perturbations of spherically symmetric BHs, previous studies suggest that, at least in specific models, the fractional deviations from GR for scalar and gravitational perturbations can be of similar orders of magnitude [110, 152], an expectation also adopted in our earlier work [94]. Comparisons between nearly static and moderately rotating BHs [153] suggest that such an expectation may remain reasonable in specific models. In particular, [152] directly compares scalar-QNM deviations with GW250114-inspired reference bands. These studies motivate our approximate comparison and help identify promising parameters for future gravitational and rotating extensions. For this purpose, we define

where each frequency is compared with its corresponding GR value at the same mass and mode indices, and the subscript “cLQBH” is added expediently to distinguish it from its GR counterpart. A positive $\delta\omega_{\mathrm{I}}$ indicates slower damping than in GR. Thus, we compare fractional deviations without equating the absolute scalar and gravitational frequencies.

Using the GW250114 constraints for the dominant $(l,m,n)=(2,2,0)$ mode [80] and translating them into the fractional deviations defined above, we adopt the rounded reference intervals $\delta\omega_{\mathrm{R}}\in[-10\%,10\%]$ and $\delta\omega_{\mathrm{I}}\in[-10\%,10\%]$. Applying the same approach to the Tian-Qin forecasts [43], we obtain the approximate ranges $\delta\omega_{\mathrm{R}}\in[-0.3\%,0.3\%]$ and $\delta\omega_{\mathrm{I}}\in[-2\%,2\%]$. In Figure 5, the light blue and light pink regions indicate deviations outside the adopted TianQin and LIGO ranges, respectively, and the pink shading takes precedence where they overlap. These shaded regions highlight parameter choices with potentially observable deviations, but do not establish observational exclusions. A more accurate assessment requires gravitational perturbations of rotating cLQBHs and a dedicated analysis of the GW data, which we leave to future work.⁵

Since the $\omega$’s in the right-hand panels of Figure 5 exhibit various types of behavior, covering both the small-deviation and the large-deviation regimes, we further list the exact values of the corresponding $\omega$ up to six decimal places in Table I. From there we confirm that the deviation from GR can indeed reach a remarkable level when the physical parameters are sufficiently large. On the other hand, by mainly referring to Figure 5, we want to carefully select another combination of $\{r_0,\mathcal{W}\}$ and extend our calculation of $\omega$ to cover the $n=1,2$ modes. Our purpose is to see the different tendencies among the different modes with a varying $\mathcal{A}$. At the same time, for this part, we will allow the deviation from GR to be sufficiently large so that the trajectory of $\omega$ in the complex plane can be clear enough (the behavior of $\omega$ for the small-deviation cases can be found in, e.g., some panels of Figure 5). Therefore, we choose $\{r_0/M,\mathcal{W}\}=\{12,2\}$ for the calculations of $\omega$ and put the results in Table II, where we have $\mathcal{A}\in[0.001,0.055]$ and we exhibit the dimensionless quantity $M\omega$ up to six decimal places.

In addition, a companion figure for this table is shown in Figure 8, exhibiting the trajectories of the three modes. Here, we show the behavior of $\omega$ for $\mathcal{A}\in[0,0.055]$, where $\mathcal{A}=0$ corresponds to the GR limit (in GR, we have $M\omega_{20}\approx0.483644-0.096759i$, $M\omega_{21}\approx0.463851-0.295604i$, and $M\omega_{22}\approx0.430544-0.508558i$ [125]). Combining the observations from Table II and Figure 8 yields a wealth of useful physical information. First of all, we notice that the changes in the $n=1$ and $n=2$ modes are more sensitive to $\mathcal{A}$ compared to the $n=0$ one, which has already been noticed in the literature [97, 154]. They can quickly and dramatically deviate from their GR counterparts in the complex plane. Secondly, it is quite clear that, for $\mathcal{A}$ in a certain regime, the $n=1$ mode

can overtake the $n=0$ mode and become the one with the longest lifetime. Evidence for this is that, in Figure 8, the red curve, which marks the mode (among the three) with the longest lifetime and the smallest damping rate, switches from the curve for $n=0$ to the $n=1$ one when $\mathcal{A}$ becomes large enough (an additional confirmation of this is shown in Figure 9). This kind of overtaking effect was noticed by some authors previously [90, 92, 155, 156]. At least for the modes considered here, the results show that the change in $\operatorname{Im}(\omega)$ is more sensitive than that in $\operatorname{Re}(\omega)$ with increasing $\mathcal{A}$. Such a tendency is consistent with the corresponding results for case (i) [84]. At the same time, we would like to emphasize that no exceptional point is observed for the parameter choices considered here, as the tracked QNM frequencies do not coalesce in the complex plane. Nevertheless, the non-Hermitian nature of the QNM problem suggests that exceptional points may (sometimes easily) be found by appropriately tuning the parameters $\{r_0, A, W\}$ [101, 157]. A systematic search for such points lies beyond the scope of this work and is left for future investigation.

$\mathcal{A}$	$M\omega_{20}$	$\delta\omega_{\mathrm{R}}$	$\delta\omega_{\mathrm{I}}$
0.001	$0.493776 - 0.100377i$	$+2.09\%$	$-3.74\%$
0.002	$0.501119 - 0.098249i$	$+3.61\%$	$-1.54\%$
0.004	$0.508980 - 0.093964i$	$+5.24\%$	$+2.89\%$
0.006	$0.513861 - 0.090959i$	$+6.25\%$	$+5.99\%$
0.008	$0.517547 - 0.088741i$	$+7.01\%$	$+8.29\%$
0.010	$0.520572 - 0.087016i$	$+7.64\%$	$+10.1\%$
0.012	$0.523168 - 0.085624i$	$+8.17\%$	$+11.5\%$
0.016	$0.527515 - 0.083496i$	$+9.07\%$	$+13.7\%$
0.020	$0.531110 - 0.081933i$	$+9.81\%$	$+15.3\%$
0.024	$0.534193 - 0.080730i$	$+10.5\%$	$+16.6\%$

Time-domain analysis

We now extend our investigation of cLQBH QNMs to the time domain by evolving the master equation (6) using the method of lines [136, 158]. We combine fourth-order spatial finite difference method (FDM) with classical fourth-order Runge-Kutta (RK4) time integration [159], and apply Sommerfeld-type boundary treatment at the boundary grid points [158]. We refer to this scheme as the FR4S method⁶. Introducing $\Pi(t,x) \equiv$

$\partial_t\Psi(t,x)$, we first rewrite the master equation as two coupled partial differential equations (PDEs),

To discretize these equations, we choose a spatial interval $[x_{\min},x_{\max}]$ and an evolution time interval $[0,t_{\max}]$, with a finite inner boundary $x_{\min}<0$ in the near-MH region and an outer boundary $x_{\max}$ in the far exterior beyond the bump. For each calculation, $t_{\max}$ is chosen sufficiently large to follow the late-time behavior of the wave function $\Psi$. We use a uniform rectangular lattice in the $(t,x)$ plane, as illustrated in [160],

encing with RK4 time integration, to reduce discretization errors and improve the accuracy of the quantitative waveform analysis below [135, 159]. where $\delta x = (x_{\max} - x_{\min})/N_x$ and $\delta t = t_{\max}/N_t$ are the spatial and time step sizes, respectively. Thus, $N_x$ and $N_t$ count the numbers of spatial and time intervals, respectively. Their values are chosen to cover the required domain and duration at the desired resolution. The indices $j$ and $k$ label the spatial grid points and time levels, respectively, which avoids confusion with the QNM overtone index $n$. At fixed $x_j$, we write $\Psi_j(t) = \Psi(t,x_j)$ and $\Pi_j(t) = \Pi(t,x_j)$ for the time-dependent fields, and $\Psi_j^k = \Psi_j(t_k)$ and $\Pi_j^k = \Pi_j(t_k)$ for their values at $t_k$. Along each grid line at fixed $x = x_j$, the potential has the time-independent value $V_j = V_{\mathrm{eff}}(r(x_j))$, where $r(x)$ is obtained numerically from Eq. (8), with a specific example shown in Figure 4.

$\mathcal{A}$	$M\omega_{20}$	$\delta\omega_R$	$\delta\omega_I$	$M\omega_{21}$	$\delta\omega_R$	$\delta\omega_I$	$M\omega_{22}$	$\delta\omega_R$	$\delta\omega_I$
0.001	$0.486596 - 0.095992i$	$+0.6\%$	$+0.8\%$	$0.402113 - 0.179523i$	$-13\%$	$+39\%$	$0.573272 - 0.229971i$	$+33\%$	$+55\%$
0.007	$0.497036 - 0.090412i$	$+2.8\%$	$+6.6\%$	$0.395466 - 0.121371i$	$-15\%$	$+59\%$	$0.603011 - 0.192020i$	$+40\%$	$+62\%$
0.013	$0.502612 - 0.086606i$	$+3.9\%$	$+10\%$	$0.394626 - 0.104457i$	$-15\%$	$+65\%$	$0.614212 - 0.182056i$	$+43\%$	$+64\%$
0.019	$0.506556 - 0.084011i$	$+4.7\%$	$+13\%$	$0.394905 - 0.094796i$	$-15\%$	$+68\%$	$0.621742 - 0.176830i$	$+44\%$	$+65\%$
0.025	$0.509673 - 0.082123i$	$+5.4\%$	$+15\%$	$0.395548 - 0.088193i$	$-15\%$	$+70\%$	$0.627542 - 0.173550i$	$+46\%$	$+66\%$
0.031	$0.512276 - 0.080685i$	$+5.9\%$	$+17\%$	$0.396323 - 0.083260i$	$-15\%$	$+72\%$	$0.632301 - 0.171299i$	$+47\%$	$+66\%$
0.037	$0.514525 - 0.079552i$	$+6.4\%$	$+18\%$	$0.397138 - 0.079370i$	$-14\%$	$+73\%$	$0.636352 - 0.169671i$	$+48\%$	$+67\%$
0.043	$0.516509 - 0.078637i$	$+6.8\%$	$+19\%$	$0.397955 - 0.076190i$	$-14\%$	$+74\%$	$0.639887 - 0.168452i$	$+49\%$	$+67\%$
0.049	$0.518288 - 0.077885i$	$+7.2\%$	$+20\%$	$0.398753 - 0.073520i$	$-14\%$	$+75\%$	$0.643025 - 0.167517i$	$+49\%$	$+67\%$
0.055	$0.519902 - 0.077257i$	$+7.5\%$	$+20\%$	$0.399526 - 0.071233i$	$-14\%$	$+76\%$	$0.645847 - 0.166788i$	$+50\%$	$+67\%$

With the lattice specified, we initialize both fields at all spatial grid points on the $k = 0$ time level using a Gaussian pulse with vanishing initial time derivative [97, 136, 161],

Here $x_c$ denotes the pulse center. We adopt $x_c/M = 5$ as our default choice in this work, which places the pulse center to the left of the bump in the tortoise coordinate, $x_c < x(r_0)$, for the parameter choices considered. To facilitate a direct comparison between the cLQBH and Schwarzschild evolutions of $\Psi$, we use the same initial data and observation prescription in both cases, thereby avoiding discrepancies arising from differences in these settings rather than from the underlying geometries [162].

We next use fourth-order spatial FDM to discretize the spatial derivative in Eq. (22). At the interior points, $2 \leq j \leq N_x - 2$, the centered stencil gives [135, 163]

The spatial truncation error is $O(\delta x^4)$. This discretization produces a system of ODEs in time. Both here and in the boundary equations below, the time derivatives are obtained by evaluating the right-hand sides of (25) through arithmetic operations on the field values at the current time level or RK4 stage.

The centered stencil is unavailable at the four outermost grid points, $x = x_0,x_1,x_{N_x-1},x_{N_x}$. To close their evolution equations, we impose the Sommerfeld relations $\partial_t\Psi= \partial_x\Psi$ at the left boundary and $\partial_t\Psi= -\partial_x\Psi$ at the right one [158]. These relations approximate purely outgoing propagation at the finite computational boundaries, thereby reducing spurious reflections. Differentiating them with respect to $t$ yields the same relations for $\Pi$. We then apply fourth-order one-sided FDM to the spatial derivatives in these Sommerfeld relations, using neighboring points toward the domain interior [163]. Writing $U_j(t) = (\Psi_j(t), \Pi_j(t))^{\mathsf{T}}$, the resulting boundary evolution equations are

Together, Eqs. (25) and (26) determine the time derivatives of both fields at every spatial grid point. At a given time level $k$, we substitute the known numerical values $\Psi_j^k$ and $\Pi_j^k$ for $\Psi_j$ and $\Pi_j$ into the corresponding right-hand sides and evaluate the time derivatives by arithmetic operations. We collect the pairs $U_j^k=(\Psi_j^k,\Pi_j^k)^{\mathsf{T}}$ for $j=0,\ldots,N_x$ vertically into a column vector $U^k$, which therefore contains both fields over the entire spatial grid at this time level. We denote the numerical operation that evaluates these right-hand sides by $\mathbf{F}$. For any supplied field-value vector $\mathbf{U}$, the output $\mathbf{F}(\mathbf{U})$ contains the corresponding time derivatives in the same order, calculated using fourth-order FDM at the interior points and Sommerfeld relations combined with fourth-order FDM at the four outermost points. This operation is applied to the current grid values and to the intermediate values constructed during RK4. Starting from the known state $U^k$, classical RK4 then advances both fields at all spatial grid points together to the next time level through [158, 163],

Here $K_1$ gives the current-state rates, while $K_2$ and $K_3$ use two estimates at the same half-step time and $K_4$ uses a full-step estimate, all constructed from $U^k$. Each evaluation of $\mathbf{F}$ uses one unchanged stage state in both the interior and boundary formulas. The final RK4 combination supplies both fields on the complete next time level, including the boundary points. Starting from $U^0$, we repeat this procedure for $k=0,\ldots,N_t-1$, recalculating all four stage-rate vectors at each step. In the end, the time-domain response will be illustrated by plots of $\ln|\Psi(t,x_{\mathrm{obs}})|$, where $x_{\mathrm{obs}}<x_{\max}$ stands for the location of the observer. For many of the calculations below, $x_{\mathrm{obs}}/M=150$ [97] provides a convenient reference location, and as long as the observer stays sufficiently far outside the bump, our major conclusions are insensitive to the precise choice of $x_{\mathrm{obs}}$.

Finally, for the nonnegative potentials considered here, the time step in the FR4S method is chosen to satisfy the following Courant-type condition for numerical stability of the interior discretization [23, 158, 163],

A typical starting choice is $\delta x/M=0.1$ and $\delta t/M=0.05$, giving $\delta t/\delta x=1/2$. This choice satisfies Eq. (28) with a safety margin in the benchmark cases and the extended parameter samples, and the condition is rechecked whenever the parameters or grid spacing change. Beyond this step-size restriction, the spatial and temporal resolutions and the boundary locations are calibrated through convergence tests and independent time-domain comparisons during the calculations. In particular, convergence can be checked by halving both $\delta x$ and $\delta t$ and comparing the resulting waveforms. These checks assess the stability and accuracy of the complete FR4S evolution, including the Sommerfeld boundary treatment, over the waveform interval of interest and also help determine the required evolution time $t_{\max}$.

We now proceed to the time-domain analysis using the FR4S method. Specifically, we shall plot $\ln|\Psi(t,x_{\mathrm{obs}})|$ [cf. (2.5)] and seek to extract physically useful information from these waveforms. Representative choices of $\{r_0,A,W\}$ will be considered to illustrate the physical features and trends, focusing mainly on $l=2$. As a first example, the parameters in the last row of Table II are chosen to allow a direct comparison with the corresponding frequency-domain results. This provides an initial link between the time-domain and frequency-domain analyses, whose complementary information can help reveal the underlying physics and yield a more complete account of the cLQBH response. For these calculations, we adopt units with $c=G=M=1$, which fixes the units of length, mass, and time in terms of the BH mass. As discussed in the previous subsection, this practical simplification leaves the dimensionless numerical results unchanged while making both the presentation and the numerical implementation more straightforward.

The temporal evolution of the wave function $\Psi$ for the cLQBH [cf. (2.5)] is exhibited in Figure 10 through $\ln|\Psi(t,x_{\mathrm{obs}})|$. Its GR counterpart is included in the same figure for comparison. As mentioned above, we adopt units with $c=G=M=1$ when solving (6) and producing this plot (so that quantities like $t$ appear to be dimensionless). This choice involves no loss of generality and simplifies both the presentation and the numerical setup in the computational codes. For Figure 10, we set $l=2$ and $\{r_0,\mathcal{A},\mathcal{W}\}=\{12,0.055,2\}$, corresponding to the last row of Table II. The observer is located at $x=x_{\mathrm{obs}}=150$ [97].

By examining Figure 10, we see that the two time-domain waveforms $\Psi(t,x_{\mathrm{obs}})$ [91], one for the cLQBH and the other for GR, almost coincide at early times, consistent with previous studies of localized potential perturbations [97]. Therefore, as shown in the figure, $t=t_{\mathrm{peak}}$ roughly marks the common maximum of the two curves. In contrast, the late-time cLQBH waveform $\Psi_{\mathrm{cLQBH}}$ exhibits a more intricate structure, with oscillation and decay patterns that differ from those of its GR counterpart $\Psi_{\mathrm{GR}}$. Here, these subscripts are introduced temporarily to distinguish the two backgrounds. More importantly, these differences are consistent with the shifts in the complex QNM frequencies $\omega$ found in the frequency-domain analysis, whose real and imaginary parts determine the oscillation periods and damping rates of the individual modes, respectively.

By adopting the so-called matrix pencil method (MPM) [165, 166, 167], we can further extract dominant QNM frequencies from the signed waveform data underlying Figure 10.⁷ For ease of reference, we denote the candidate complex frequencies by $\hat{\omega}_{lk}$. The index $k$ distinguishes their labeling from the frequency-domain branch index $n$, so no correspondence between them is assumed before comparison. After obtaining the candidate frequencies, we vary the fitting settings to screen out spurious roots and retain stable candidates, using a consistency check analogous to those employed in the matrix method [128] and the CFM [169]. Following this screening, the two dominant modes for Figure 10 are found to be $\hat{\omega}_{20}\simeq0.399526-0.071233i$ and $\hat{\omega}_{21}\simeq0.519903-0.077257i$. A much weaker component consistent with the frequency-domain $n=2$ mode is also recovered, but its frequency is not quoted because it is extracted with lower precision. Here the retained modes with positive real parts are listed in order of increasing $\lvert\mathop{\mathrm{Im}}\hat{\omega}_{lk}\rvert$, so a smaller $k$ corresponds to a longer lifetime, with $k=0$ labeling the longest-lived mode [46]. Comparison with Table II identifies $\hat{\omega}_{20}$ with $\omega_{21}$ and $\hat{\omega}_{21}$ with $\omega_{20}$, and the differences in both the real and imaginary parts are smaller than $10^{-5}$. Whenever small differences arise, the frequency-domain values serve as the reference because they have undergone systematic convergence tests and independent cross-checks. Such differences can result from the finite fitting window and non-QNM contributions in waveform-based extraction, particularly for weak components [167, 170]. The prominence of the longer-lived modes is not surprising, since the late-time fitting window favors the less damped excited QNMs, whose contributions decay more slowly during ringdown while the more strongly damped components fade. In particular, the continued $n=1$ branch outlives the $n=0$ branch, providing further time-domain evidence of the overtaking already identified in the frequency-domain analysis [90].

As a complementary check of the frequencies identified among the candidates $\hat{\omega}_{lk}$ against the waveform underlying Figure 10, we reconstruct the waveform using the modes retained after numerical screening [171]. The initial MPM model contains $p$ exponential components. After screening out spurious candidates by varying the fitting settings, each retained positive-real-frequency mode $\hat{\omega}_{lk}$ is considered together with its companion $-\hat{\omega}_{lk}^{*}$. The positive-real-frequency members of these pairs are numbered $k=0,\ldots,q$ in order of increasing $\lvert\mathop{\mathrm{Im}}\hat{\omega}_{lk}\rvert$, in accordance with the convention above. Thus $q$ is the largest index in the retained set, which contains $q+1$ mode pairs, so that $2(q+1)\leq p$. For the oscillatory modes considered here, the general form of the reconstructed real-valued waveform is [172]

where the superscript $*$ denotes complex conjugation. Each complex coefficient $C_{lk}$ determines the amplitude and phase of the corresponding mode at $t=t_a$, the beginning of the fitting window [167]. The retained modes need not have comparable extraction accuracy or amplitudes. In practice, the reconstruction may therefore use only those modes determined with sufficient accuracy and carrying appreciable contributions, without necessarily including every term in Eq. (29). In the late-time fitting window considered here, the less damped modes also have larger relative amplitudes and are more accurately extracted, so these considerations favor the same dominant modes.

To determine these coefficients, we fit the reconstructed waveform in Eq. (29) to the original FR4S data at the sampling times $t=t_a+j\delta t_s$. For this purpose, the sampled values and the unknown coefficients are collected into the column vectors

where T denotes transposition. At each sampling time, Eq. (29) expresses the reconstructed value as a linear combination of the coefficients in C. The exponential factors multiplying these coefficients therefore form a numerical matrix E, whose elements are

where $j=0,\ldots,N_s$ and $k=0,\ldots,q$. Consequently, EC gives the reconstructed values at the sampling times, and the fitting relation becomes

The matrix E and the data vector y consist entirely of known numbers, obtained from the retained MPM frequencies and the original FR4S samples. With the frequencies $\hat{\omega}_{lk}$ held fixed, we determine the coefficients $C_{lk}$ by minimizing

Here $(\mathrm{EC})_{j+1}$ is the reconstructed waveform value at $t=t_a+j\delta t_s$, and $y_j$ is the corresponding FR4S value defined in Eq. (30). This is a linear least-squares fit [167, 170]. Only the coefficients are fitted, and each $C_{lk}^{*}$ is constrained to be the complex conjugate of $C_{lk}$. For plotting, as mentioned above, we may further retain only the terms with sufficiently accurate frequencies and appreciable amplitudes in Eq. (29), keeping their fitted coefficients unchanged. This truncation gives an approximate reconstruction, without implying that the omitted modes are absent. Its accuracy is assessed by comparing the resulting waveform directly with the FR4S data. We note that this coefficient fit, together with the other linear algebra operations in MPM, can be carried out using the built-in functions of, e.g., Mathematica.

Substituting the resulting numerical coefficients into Eq. (29), with the sum restricted to the included

modes when appropriate, gives an explicit reconstructed waveform as a function of $t$. This waveform can be compared directly with the original FR4S waveform, without an additional normalization. Accordingly, for the case shown in Figure 10, we plot in Figure 11 the FR4S waveform $\Psi(t,x_{\mathrm{obs}})$ together with its reconstruction from Eq. (29), where $t$ is expressed in dimensionless form through the units with $c=G=M=1$ adopted earlier. The two curves agree very well, even though only two modes, $\hat{\omega}_{20}$ and $\hat{\omega}_{21}$, are used for the reconstruction. In fact, a closer examination of the contributions of the different modes, measured by $|C_{2k}|/\max\{|C_{20}|,|C_{21}|,\ldots,|C_{2q}|\}$, shows that these two modes make comparable contributions to $\Psi_{\mathrm{rec}}$, while each remaining fitted component contributes less than $0.03\%$ of the largest one. This explains why these two modes are sufficient to provide an adequate truncated reconstruction over the interval shown. On the other hand, the comparable contributions of these two modes suggest that their interference explains the more intricate oscillation and decay patterns in Figure 10, consistent with similar effects observed in other multi-mode ringdown waveforms [168].

The preceding example illustrates a waveform in which two modes make comparable contributions, giving the late-time behavior of $\Psi$ a more intricate structure (cf. Figure 10). We now examine how the bump amplitude and width affect the response at a fixed location, using the four configurations shown in Figure 12 and Figure 13. We choose $r_0=10$, an intermediate location between the nearer and farther bumps considered above, and compare four combinations of the amplitude $\mathcal{A}$ and the width $\mathcal{W}$. The four resulting bump profiles, shown in the insets of Figure 12, combine $\mathcal{A}=0.003$ and $0.03$ with $\mathcal{W}=0.5$ and $2.5$. The calculations retain $l=2$, $x_{\mathrm{obs}}=150$, the same initial data, and the units $c=G=M=1$. All four MPM extractions use the common fitting window $[t_a,t_b]=[230,290]$ and retain the previous settings $\delta t_s=0.2$, $p=10$, and $l=100$. Figure 13 compares the corresponding FR4S and reconstructed waveforms over $t\in[225,295]$, with the vertical dashed lines marking the fitting-window boundaries.

As seen in Figure 12, the smallest and narrowest bump produces only a slight departure from GR, qualitatively consistent with the small frequency shifts at low amplitudes shown in Figure 5. Increasing either the amplitude or the width reduces $\operatorname{Re}\omega_{20}$ across these four configurations, giving a longer oscillation period. The damping behavior is less regular. Three of the four configurations decay faster than in GR, while the fourth decays more slowly, as reflected by their oscillation envelopes lying below or above the GR envelope. In particular, increasing the width enhances the damping at the smaller amplitude but reduces it at the larger amplitude. The largest departure occurs for the larger and broader bump, for which $\operatorname{Re}\omega_{20}$ and $|\operatorname{Im}\omega_{20}|$ are approximately $6.9\%$ and $23.8\%$ smaller than their GR counterparts, respectively. This again confirms the finding in Sec. III A that the change in $\operatorname{Im}\omega$ is more sensitive than that in $\operatorname{Re}\omega$. These results show the combined influence of the amplitude and width, without establishing that either alone controls the response. Figure 13 shows that all four late-time waveforms are accurately reconstructed using a single dominant mode, consistent with the usual dominance of the fundamental mode during ringdown [46], since the more strongly damped overtones can become much weaker at sufficiently late times [173]. Frequency-domain continuation identifies this dominant mode with the continued $n=0$ branch, i.e., the $\omega_{20}$ mode, and the MPM and shooting results agree when rounded to six decimal places. Thus, a noticeable departure from GR does not necessarily involve competition between two comparable mode contributions, although more distant bumps can allow other branches to overtake the continued $n=0$ mode [90]. Moreover, the agreement in Figure 13 extends slightly beyond the fitting window without readjusting the coefficients. Together with the frequency-domain verification, this supports the physical interpretation of the reconstruction in terms of QNM ringing and motivates further investigation of the bump location and the timescale on which the waveform differences become appreciable.

Time for deviation and the observational window

The waveform comparisons in Figure 12 illustrate how the bump amplitude and width affect the time-domain response at a fixed $r_0$. We now turn to the bump location $r_0$, which changes the separation between the localized deformation (viz., the bump) and the main potential barrier and has received particular attention in studies of the bump effect in recent years [90, 91, 97]. These studies show that moving a small bump outward to the positive-$x$ direction can substantially modify the QNM spectrum while leaving the early ringdown nearly unchanged (cf., Figure 10 and Figure 12), motivating a complementary examination of the waveforms. In our cLQBH model, the spiral trajectory in Figure 6 further illustrates the nontrivial spectral dependence on $r_0$ and raises the question of how this dependence manifests itself in the time-domain response. We therefore extend the preceding analysis by holding $A$ and $W$ fixed and varying $r_0$. By studying the bump effect in cLQBHs, we connect the general question of how localized potential deformations affect BH ringing with the dynamical consequences of quantum-corrected geometries. The dependence on $r_0$ provides a natural point of contact, allowing us to examine how the location of the correction shapes the time-domain response. Accordingly, the choices of $\{A,W\}$ and the ranges of $r_0/M$ warrant careful consideration, both to remain within the physically admissible regime and to reveal waveform trends that can be clearly identified and interpreted. These considerations guide the following waveform comparisons using the FR4S method described in Sec. III. To give a set of examples, we calculate $\Psi(t,x_{\rm obs})$ for three representative values of $r_0$, keeping $\{A,W\}=\{0.004,1.25\}$ fixed. The resulting waveforms are shown

in Figure 14, where we adopt $l=2$, $x_{\rm obs}=150$, and units with $c=G=M=1$, as before. These comparisons show that varying the bump location changes how the cLQBH waveform departs from its GR counterpart, complementing the nontrivial dependence of the QNM spectrum on $r_0$ illustrated in Figure 6. Specially, at least for the late-time windows, we notice again that the departure from that of GR is more explicit with a relatively larger $r_0$.

More importantly, these waveforms draw attention to a quantity of particular interest, the approximate time at which the difference between the two responses becomes appreciable, which we denote by $t_{\rm dev}$. The onset of this appreciable departure can be estimated visually, and similar behavior is evident in previous waveform comparisons [91, 97]. These studies also relate the delayed response to the location of the potential perturbation through wave propagation and scattering. Building on these results, a systematic examination of $t_{\mathrm{dev}}$ in the present cLQBH model can further clarify how the bump location controls the emergence of waveform differences. The examples in Figure 14 suggest that $t_{\mathrm{dev}}$ need not vary monotonically with $r_0$, motivating a more detailed examination of the full dependence. We therefore undertake a semi-quantitative study of this behavior. Although the departure itself is a physical feature of the waveforms, assigning a reproducible time to its onset requires an explicit criterion. Such a criterion allows us to identify $t_{\mathrm{dev}}$ automatically and consistently across different parameter choices, while retaining its connection to visual inspection. For this purpose, we combine the normalized waveform distance used in [174] with Gaussian weighting [175] to introduce the following local measure of waveform deviation,

where $t_k=k\delta t$, $j$ labels the time steps relative to $t_k$, and $T$ sets the width of the Gaussian weight. This dimensionless quantity measures the local difference between the two waveforms around any given $t_k$, with the denominator providing a common normalization. For a nonzero denominator, this measure satisfies $0\leq D_T\leq\sqrt{2}$. We evaluate the sums directly from the uniformly sampled FR4S data and examine $D_T(t_k)$ chronologically, starting at the common GR amplitude peak $t_{\mathrm{peak}}$ (indicated in Figure 14). For the cases considered here, $D_T$ is initially below a prescribed threshold $\delta$. We locate the first sample satisfying $D_T(t_k)\geq\delta$, whose preceding sample therefore satisfies $D_T(t_{k-1})<\delta$, and determine the crossing time by a straightforward linear interpolation,

Thus, interpolation is applied only to the two adjacent diagnostic values that bracket the threshold, while the waveforms entering Eq. (34) are used directly. The corresponding continuous expression replaces the sums in Eq. (34) with integrals over $[t_k-6T,t_k+6T]$. Accordingly, we set $J=6T/\delta t$, so that the discrete and continuous expressions cover the same interval. This truncation retains the Gaussian weight out to six standard deviations on either side of its center, beyond which the weight is negligible. Checks extending the interval to $\pm8T$ confirm a negligible change in the extracted deviation times for the tested cases. Within this interval, we choose $T=5M$, a width parameter comparable to half the reference GR ringing period, approximately $6.5M$. Gaussian widths on this half-period scale have also been considered in local oscillation analysis [176], providing a methodological precedent for this choice. With the sampling interval $\delta t=0.05M$ [cf., the context around (28)], we therefore have $J=600$.

We adopt the fixed threshold $\delta=0.075$, which is of the same order as the value $0.1$ used for a similarly normalized distance in electrical waveform comparisons [174]. This provides a reference scale rather than a prescribed threshold for the present problem. We assess the suitability of our choices by examining the sensitivity of the extracted times to nearby values of $T$ and $\delta$, checking numerical resolution and the summation interval, and comparing the resulting markers with the waveforms themselves. These checks support the use of the adopted criterion for our semi-quantitative analysis, while the precise onset times remain dependent on the stated criterion. We keep the same $T$ and $\delta$ for all bump locations to ensure a consistent comparison.

Applying this procedure to Figure 14 gives the deviation times marked by the vertical steel-blue dashed lines, with their numerical values indicated alongside them. These examples again show that $t_{\mathrm{dev}}$ does not vary monotonically with $r_0$. To illustrate the time interval entering Eq. (34), the pale-yellow shading extends over $t_{\mathrm{dev}}\pm6T$, with its intensity proportional to the Gaussian weight. Its gradual fading shows how samples farther from the central time are assigned smaller weights. As a more direct visual cross check, Figure 15 presents a local enlargement of the signed waveforms for the $r_0/M=33$ case in the bottom panel of Figure 14, with the axis ranges and aspect ratio adjusted appropriately to make their separation easier to inspect. The marker at $t_{\mathrm{dev}}/M\approx218.2$ lies near the transition from nearly overlapping waveforms to an increasingly visible separation. This agreement supports interpreting the extracted time as an approximate onset of appreciable deviation, consistent with the local averaging used in the definition of $t_{\mathrm{dev}}$.

Beyond their phenomenological significance, the results discussed above lead us to consider their observational implications. In this work, we aim not only to connect the cLQBH model with the bump effect on QNMs, but also to explore whether the resulting phenomena could have observable consequences. Existing studies have examined several complementary aspects of QNM observability, including whether the ringdown frequencies fall within the sensitive frequency bands of planned detectors [43] and whether multiple modes can be detected and resolved for black-hole spectroscopy [83, 177], etc.

Within the present framework, there are two closely related ways to investigate this question. The first is to calculate the QNM frequencies $\omega$ and compare them with their GR counterparts through, for example, the fractional deviations defined in Eq. (21). Comparing these deviations with the measurement accuracies of current and planned detectors, including LVK, LISA, and TianQin, allows one to assess possible constraints on the parameters controlling the QNMs of a gravitational theory or model. Such constraints can be investigated using existing observations or forecasts in which GR is assumed to remain consistent with future data [43, 80]. Our comparisons in Figure 5 and Table I provide an exploratory assessment in this direction. Related studies have also considered this approach using scalar QNMs [94, 152]. The second is to investigate whether the bump-induced changes in the QNMs and waveforms become sufficiently pronounced while the signal still retains an appreciable amplitude. In particular, even substantial changes in the QNM spectrum may leave the early ringdown nearly unchanged, with more apparent waveform differences emerging only later [91, 97]. The observational relevance of the bump effect has already been discussed, with previous studies emphasizing that small corrections to the early ringdown and weak late-time signals can limit the observational accessibility of the modified QNM spectrum [91, 106]. These considerations motivate a more systematic examination of when waveform deviations emerge relative to signal decay in the present cLQBH model. Although closely connected to the first approach, this second approach explicitly considers whether the signal persists with sufficient amplitude for the deviations to become apparent. In this section, we pursue this direction by comparing the deviation time with the signal decay, thereby introducing an operational observation window for the scalar waveforms calculated above.

In practice, assessing the observational significance of BH ringing requires selecting a suitable time interval for the signal analysis [72, 178]. A natural and practical starting point is the amplitude peak, where the signal is strongest, so that the interval includes the strong part of the response while excluding the initially weak signal before its main arrival. The time of such a peak also serves as a reference time in gravitational ringdown studies, including analyses that incorporate overtones [73, 164, 173]. Accordingly, we adopt the time of the common GR reference peak, $t=t_{\rm peak}$, introduced above, as the starting point of our waveform comparison, since the cLQBH waveform reaches its peak at approximately the same time (see, e.g., Figure 10). For example, for $\{r_0/M,\mathcal{A},\mathcal{W}\}=\{18,0.004,1.25\}$ in the middle panel of Figure 14, the estimated peak-time difference is $\lvert t_{\rm peak,cLQBH}-t_{\rm peak,GR}\rvert/M\simeq0.0004$. After the peak, the signal generally decays, and deviations that emerge sufficiently late may therefore occur when its amplitude has already become very small. This motivates introducing an approximate faintness threshold to determine the endpoint of the interval, in the spirit of the relative-amplitude cutoffs used in numerical ringdown analyses [164]. We denote the first post-peak time at which a suitable amplitude measure falls below this threshold by $t_{\rm faint}$. The interval between $t_{\rm peak}$ and $t_{\rm faint}$ then provides a simple framework for examining whether the bump-induced deviations emerge while the signal retains an appreciable amplitude. Although a quantitative assessment of GW detectability requires gravitational perturbations, source properties, and the detector response and noise, the scalar waveforms obtained here allow us to explore the interplay between waveform deviation and signal decay and its dependence on the cLQBH parameters. Beyond providing physical insight into the bump effect, this analysis can serve as a prototype for more complete investigations of its observational implications.

Relative-amplitude cutoffs have been used to determine the endpoint of numerical ringdown fits. In particular, Berti et al. truncate each multipolar waveform when its amplitude falls below $10^{-3}$ or $10^{-4}$ of its own peak amplitude to limit late-time numerical contamination [164, 167]. Following this relative-amplitude prescription, we introduce a fixed fraction $\eta$ of the approximate peak amplitude of each cLQBH waveform. Writing $\Psi_{\mathrm{cLQBH}}$ for the waveform denoted by $\Psi$ above, we define $\Psi_{\mathrm{peak}} \equiv\lvert\Psi_{\mathrm{cLQBH}}(t_{\mathrm{peak}},x_{\mathrm{obs}})\rvert$, using the common reference time $t_{\mathrm{peak}}$, which approximately marks the peaks of both the cLQBH and GR waveforms in Figure 10. For each parameter choice, $t_{\mathrm{faint}}$ is the first post-peak time at which the local amplitude measure defined below falls below $\eta\Psi_{\mathrm{peak}}$. The chosen $\eta$ provides a reasonable and practical reference for comparing the waveforms, rather than an absolute boundary beyond which the signal becomes unobservable for any realistic analyses.

To determine $t_{\mathrm{faint}}$ systematically, we estimate the local amplitude by averaging the absolute waveform values at two successive stationary points. Here we use $t_k$ to label successive stationary points (local maxima and minima identified approximately from the sampled waveform) of $\Psi_{\mathrm{cLQBH}}(t,x_{\mathrm{obs}})$, ordered chronologically. These selected times are generally not equally spaced. We define

The average is assigned to the earlier extremum time $t_k$. Using the extrema avoids interpreting individual zero crossings of the oscillatory waveform as a loss of signal amplitude.

As illustrated in Figure 10, the cLQBH and GR waveforms nearly coincide around their maxima, so the fixed reference time $t_{\mathrm{peak}}$ provides a good approximation to their common peak time. Starting from this reference, we identify the first extremum time $t_k$ at which $\overline{\Psi}(t_k) \le\eta\Psi_{\mathrm{peak}}$. For the cases considered here, the initial post-peak estimate exceeds the threshold, so the preceding estimate satisfies $\overline{\Psi}(t_{k-1}) > \eta\Psi_{\mathrm{peak}}$. We then determine $t_{\mathrm{faint}}$ by linear interpolation between these two amplitude estimates,

which is similar to the treatment of (35). This gives the first threshold crossing of the linearly interpolated amplitude measure, allowing $t_{\mathrm{faint}}$ to lie between extrema rather than restricting it to a stationary point. We evaluate this condition separately for each cLQBH waveform, using its own approximate peak amplitude $\Psi_{\mathrm{peak}}$. If the condition is never satisfied within the available data, no $t_{\mathrm{faint}}$ is reported for that evolution.

We adopt the representative value $\eta= 0.003$, corresponding to an average amplitude of 0.3% of the approximate waveform peak. This is of the same order as the $10^{-3}$ relative-amplitude cutoffs discussed above and provides a practical reference for identifying a substantially decayed signal. Keeping $\eta$ fixed allows a consistent comparison across the cLQBH parameter choices. As with $t_{\mathrm{dev}}$, the extracted time depends on the stated criterion, and its sensitivity to $\eta$ can be assessed by varying the threshold while keeping the waveform data and amplitude prescription fixed. The interval $t \in[t_{\mathrm{peak}},t_{\mathrm{faint}}]$ then defines an “observation window” anticipated in Sec. I. We use this window to assess whether appreciable departures from GR emerge while the cLQBH signal retains an appreciable relative amplitude. In particular, $t_{\mathrm{peak}} \le t_{\mathrm{dev}} < t_{\mathrm{faint}}$ indicates that the deviation criterion is met before the adopted faintness threshold is reached. We now compare $t_{\mathrm{dev}}$ and $t_{\mathrm{faint}}$ as the bump parameters vary.

In Figure 16, we present our first example illustrating the behavior of $t_{\mathrm{dev}}$ for the case discussed in Figure 14, with the curves of $t_{\mathrm{faint}}$ and $t_{\mathrm{peak}}$ also shown for reference. The light-blue shaded region marks the time interval between $t_{\mathrm{peak}}$ and $t_{\mathrm{faint}}$. According to the definitions and discussion given above, a $t_{\mathrm{dev}}$ lying within this light-blue region indicates that the associated bump effect is more likely to be observable from the QNM point of view. It is therefore clear that a “bowl-shaped” region, enclosed by the curves of $t_{\mathrm{dev}}$ and $t_{\mathrm{faint}}$, exists between $r_0 \approx10.6$ and $r_0 \approx21.3$ (in units where $c=G=M=1$). Within this region, the deviation between the cLQBH and GR signals occurs before the signal becomes too faint to be noticeable (based on our current definitions). This result provides a clue to, and an explanation of, the fact that a large deviation of the QNM from its GR counterpart alone does not guarantee its detectability. Instead, owing to the late-time features of the bump effect, we must also ask whether the corresponding $t_{\mathrm{dev}}$ can enter the favored region of possible detection, as illustrated by the light-blue region in Figure 16. We therefore use this criterion to define the observational window mentioned above. Specifically, the bump effect can be employed to search for geometries whose parameters yield a $t_{\mathrm{dev}}$ curve that partially enters the light-blue shaded region shown in Figure 16. Otherwise, it falls outside the observational window, even if the resultant QNM exhibits a sufficiently large deviation from its GR counterpart and is thus distinguishable by detectors at an adequate signal-to-noise ratio.

As expected, varying $\mathcal{A}$ and $\mathcal{W}$ can change the extent of the bowl-shaped region discussed above. In particular, for certain combinations of $\{\mathcal{A},\mathcal{W}\}$, the $t_{\mathrm{dev}}$ curve may lie entirely above the $t_{\mathrm{faint}}$ curve, indicating that the waveform deviations become appreciable only after the signal has fallen below the adopted faintness threshold. An example is shown in the upper panel of Figure 17. Conversely, suitable parameter choices can enlarge the bowl-shaped region, allowing the bump effect to become potentially observable over a wider range of $r_0$, as illustrated in the lower panel of Figure 17. Since localized potential bumps may also arise from matter surrounding a BH [90], these results may have implications for environmental effects on BH ringing. In particular, the location of such matter could play an important role in whether the associated waveform deviations become appreciable before the signal becomes too faint.

Building on the above results, we now address several key points and discuss their implication in turn. First, the observation window identified above does not imply that every parameter choice within the bowl-shaped region produces QNM frequency shifts large enough to be resolved by current or future detectors. As discussed earlier in this section, the comparison of $t_{\mathrm{dev}}$ with $t_{\mathrm{faint}}$ and the fractional frequency deviations $\delta\omega_{R,I}$ in (21) provide complementary assessments of observational relevance. More generally, a detectable departure from GR requires the relevant frequencies to fall within the detector’s sensitive band, the waveform differences to be dis-

tinguishable at the available signal-to-noise ratio after accounting for parameter degeneracies, and the deviations to emerge while sufficient signal remains [47, 83, 177]. Favorable conditions in only some of these respects need not ensure detectability. Accordingly, we assess the spectral and temporal criteria jointly through representative examples, leaving systematic parameter scans, consistent gravitational-perturbation calculations, and dedicated GW data analysis to future work. For example, the parameter choice $\{r_0, A, \mathcal{W}\} = \{18, 0.004, 1.25\}$ within the bowl-shaped region of Figure 16 yields a dominant scalar QNM frequency $M\omega_{20} \simeq0.495349 - 0.088255\mathrm{i}$, corresponding to $\delta\omega_R \approx+2.42\%$ and $\delta\omega_I \approx+8.79\%$ according to (21). Both deviation magnitudes also exceed the approximate TianQin reference accuracies adopted around Figure 5. This example demonstrates that appreciable frequency shifts can coexist with waveform deviations emerging within the adopted observation window, providing a concrete motivation for considering the spectral and temporal criteria together when assessing the observational relevance of a QNM problem. It is therefore important to consider both aspects when evaluating the detectability of a QNM signal.

This brings us to a more careful examination of the combined assessment of the detectability of the QNM frequency shifts in our cLQBH model, although a more complete investigation of this kind should be left to future work, as mentioned above. As an illustrative example, the parameter choice $\{r_0/M, A, \mathcal{W}\} = \{20, 0.0012, 3\}$ leads to $M\omega_{20} \simeq0.477956 - 0.084511\mathrm{i}$. According to (21), the corresponding fractional deviations from GR are $\delta\omega_R \approx-1.2\%$ and $\delta\omega_I \approx+12.6\%$. Using the approximate TianQin reference accuracies introduced around Figure 5, both deviations exceed the adopted reference intervals and are therefore potentially resolvable by TianQin, subject to the observational qualifications discussed above. In contrast, for case (i), applying the same criteria and assigning the same amplitude to $\delta V_i$ yields only $\delta\omega_R \approx+0.02\%$ and $\delta\omega_I \approx-1.97\%$. Both deviations are smaller in magnitude than the adopted TianQin reference accuracies and would therefore remain unresolved by this criterion. This specific example illustrates a crucial point: by taking advantage of the bump effect, one can potentially amplify the spectral signatures of quantum corrections to LQBHs and their QNMs through the enhanced response of the spectrum to localized potential deformations, and importantly, despite the special form of $\chi(r)$ adopted here, this amplification is not peculiar to that profile but reflects a generic feature of the bump effect, as argued in previous studies [90, 91, 97]. This possibility is particularly relevant because, in many LQBH models, the quantum corrections admitted by the underlying theory or arising in physically natural parameter regimes produce only weak QNM deviations from GR, making these deviations difficult to identify observationally [94, 110]. Although the holonomy profile in case (ii) is phenomenological and is primarily intended to produce a simple localized potential deformation whose location, amplitude, and width can be controlled separately, it is nevertheless an admissible choice within the background framework considered here. As pointed out in the preceding section, this makes it a suitable first step in this direction. More importantly, it demonstrates a possible and potentially viable route for amplifying the QNM signatures of quantum corrections to BH geometries, thereby opening an accessible observational window for distinguishing the relevant QNMs from their GR counterparts. Further applications of this mechanism merit deeper investigation.

To further investigate the generic character of the mechanism illustrated in Figure 16 and Figure 17, we consider an illustrative odd-parity QNM calculation using the metric-based effective potential employed in Refs. [108, 110] and in the dressed-metric treatment of Ref. [116]. Here, we incorporate quantum corrections through the background metric but do not include possible additional corrections to the perturbation dynamics. This provides an effective odd-parity model whose consistency with the complete cLQBH gravitational dynamics remains to be established. Within this prescription, the master equation retains the form of (6), with $\Psi$ replaced by $\Psi_{\mathrm{odd}}$, $V_{\mathrm{eff}}$ replaced by $V_{\mathrm{odd}} = V_{\mathrm{RW}} + \delta V_{\mathrm{odd}}$, and the tortoise coordinate still defined by (8). The two contributions to $V_{\mathrm{odd}}$ are given by

Here, $V_{\mathrm{RW}}$ is the Schwarzschild Regge–Wheeler (RW) potential [110], and $\delta V_{\mathrm{odd}}$ plays the role of the localized correction in this odd-parity problem. For the configurations considered below, it contains a negative trough followed by a positive peak, reversing the ordering found in the scalar bump. Following the conventions in Figure 1, we characterize this correction by its zero-crossing location $\hat{r}_0$, its dimensionless peak-to-trough amplitude $\hat{A}$, and its dimensionless peak-to-trough separation $\hat{\mathcal{W}}$. These quantities are determined from $\delta V_{\mathrm{odd}}$ and generally differ from the corresponding scalar quantities for the same background, although, as will be seen shortly, they can be quite close numerically. We choose $\{\hat{r}_0,\hat{A},\hat{\mathcal{W}}\} = \{16, 0.006, 1.3\}$ as a representative example, adopting $c = G = M = 1$. The corresponding scalar parameterization of the same geometry is $\{r_0,A,\mathcal{W}\} \approx\{16.11, 0.00602, 1.298\}$. Thus, the two parameterizations are numerically close in this example. The top panel of Figure 18 compares $V_{\mathrm{odd}}$ with $V_{\mathrm{RW}}$ for $l = 2$, while the middle panel shows the corresponding FRS waveforms at $x_{\mathrm{obs}} = 150$, using the same initial data and numerical procedures as for the scalar-perturbation case. The middle panel shows initially close agreement with the GR waveform, denoted by $\Psi_{\mathrm{RW}}$, followed by an appreciable deviation at later times, consistent with previous time-domain studies of a RW potential supplemented by a localized bump [91]. and with our scalar results (see, e.g., Figure 10). The extracted complex QNM frequency for the $\{l,n\}=\{2,0\}$ mode is $M\omega_{20}^{\mathrm{odd}}\simeq0.363764-0.081773i$. The superscript “odd” distinguishes this frequency from the scalar results. Comparing this result with the Schwarzschild value $M\omega_{20}^{\mathrm{RW}}\simeq0.3736716844-0.0889623157i$ [179] using (21) yields $\delta\omega_{\mathrm{R}}\approx-2.65\%$ and $\delta\omega_{\mathrm{I}}\approx+8.08\%$. Both magnitudes comfortably exceed the approximate TianQin reference accuracies adopted around Figure 5. This result helps motivate future investigations with consistent gravitational waveforms as a promising extension of the present analysis. More importantly, to examine the temporal criterion, we define $\hat{t}_{\mathrm{dev}}$ and $\hat{t}_{\mathrm{faint}}$ using the same prescriptions as (35) and (37), now comparing $\Psi^{\mathrm{odd}}$ with $\Psi^{\mathrm{RW}}$ and determining the faintness time separately from each odd-parity waveform. The reference time $\hat{t}_{\mathrm{peak}}$ is obtained from the RW waveform. For the example above, $\hat{t}_{\mathrm{peak}}\approx146.7$ and $\hat{t}_{\mathrm{dev}}\approx202.6$, as indicated in the middle panel. Keeping $\hat{A}=0.006$ and $\hat{W}=1.3$ fixed while varying the actual zero-crossing location $\hat{r}_0$ produces the bottom panel. A bowl-shaped region again appears, with $\hat{t}_{\mathrm{dev}}<\hat{t}_{\mathrm{faint}}$ between the intersections at $\hat{r}_0\approx10.2$ and $\hat{r}_0\approx24.0$. The representative case $\hat{r}_0=16$ therefore satisfies both the temporal and spectral criteria adopted here. The persistence of this behavior in the effective odd-parity calculation supports its relevance beyond the scalar examples. Although it does not establish universality for arbitrary backgrounds or potential deformations, this similarity between two distinct perturbation problems suggests that analyses based on the deviation time offer a promising direction for future detectability studies.

Finally, the results presented here highlight the bridge that our work builds between the investigation of QNMs in LQBHs and the generic bump effect. On the one hand, as discussed above, our analysis serves as a prototype for a more complete framework and provides preliminary evidence that the bump effect can aid the prospective extraction of LQBH QNMs by offering a potentially viable route to amplify the spectral signatures of quantum corrections to the geometry and the associated QNM deviations from GR, thereby increasing the prospect that these signatures fall within the reach of the detectors considered here. Moreover, Figure 16, Figure 17, and Figure 18 suggest a selection mechanism. For a localized correction to the background geometry, whether induced by an effective geometric deformation or by surrounding matter [90, 106, 107], observable QNM signatures may be favored only within certain parameter ranges, even when other parameter choices produce larger departures from GR. This behavior reflects a competition between the magnitude of the departure from GR and the time at which the waveform deviation emerges. Under the criteria adopted here, parameter choices that achieve a favorable balance between these two effects are more likely to yield observationally relevant QNM signals. On the other hand, our case (ii) [cf. (1.1)] provides a concrete BH platform for investigating the bump effect. Several characteristic features identified in previous studies of prescribed potential deformations are recovered in this setting, including spiral trajectories of the QNM frequencies in the complex plane, mode overtaking, and waveforms in the time domain that initially remain close to their GR counterparts before departing at later times [90, 91, 97, 149], as illustrated by Figure 6 and Figure 12, together with Table II. The cLQBH results further illustrate the value of combining spectral and temporal criteria in assessing the observational relevance of the bump effect, providing a preliminary step from theoretical modeling toward observational applications. The flexibility of case (ii) also allows other aspects of the bump response to be examined, including its amplitude, width, and location dependences, as illustrated in Figure 5, Figure 7, and Figure 16.

Conclusions

In this paper, we have investigated the QNM spectrum and time-domain response of a massless, minimally coupled test scalar field on static, spherically symmetric cLQBH backgrounds. Our main focus has been a phenomenological holonomy function with a Gaussian-localized radial profile, with the previously studied decreasing profile serving as a reference. The localized holonomy correction modifies both the scalar effective potential and the relation between the areal radius and the tortoise coordinate. Consequently, this model provides a realization of the bump effect in which the scattering potential and wave-propagation geometry follow consistently from the same quantum-corrected background.

The resulting potential deformation contains a positive peak and a negative trough. We have characterized it through its dimensionless location $r_0/M$, dimensionless amplitude $\mathcal{A}$, and dimensionless width $\mathcal{W}$, allowing the effects of these quantities to be examined separately. Our calculations concentrate on configurations for which the total effective potential remains positive throughout the black hole exterior. Within this parameter regime, localized corrections produce a variety of spectral behaviors. Varying the bump location generates spiral trajectories of the QNM frequencies in the complex plane, while changes in the amplitude and width affect both the oscillation frequencies and damping rates. The response depends on the combined parameter choices, and the damping can either increase or decrease relative to Schwarzschild.

Tracking the $l=2$ modes continuously from their Schwarzschild limits reveals that the branches associated with the first two overtones are substantially more sensitive to the localized correction than the continued fundamental branch. In particular, the continued $n=1$ branch can overtake the continued $n=0$ branch and become the longest-lived among the three tracked modes. Here, the overtone indices identify continuation from Schwarzschild rather than an ordering by the damping rate at each parameter value. This change in lifetime ordering occurs without a coalescence of the tracked complex frequencies and no exceptional point was identified in the cases examined. Our time-domain calculations complement these spectral results. For representative configurations, the early cLQBH response remains close to its Schwarzschild counterpart, whereas the subsequent ringing exhibits different oscillation and decay patterns. Frequencies extracted using the matrix pencil method agree closely with the frequency-domain shooting results, and reconstructions using the dominant extracted modes reproduce the numerical waveforms over the selected ringdown intervals. These comparisons also confirm the mode-overtaking behavior: in the representative two-mode example, the longer-lived contribution corresponds to the branch continued from the first Schwarzschild overtone. Other configurations are accurately described by a single dominant mode, demonstrating that appreciable waveform changes do not necessarily require competition between comparably excited modes.

To examine when these differences become significant relative to signal decay, we have introduced a deviation time $t_{\mathrm{dev}}$, defined through a Gaussian-weighted, normalized waveform distance. We compare it with a faintness time $t_{\mathrm{faint}}$, determined by a prescribed relative-amplitude threshold. Together with the reference peak time, these quantities define an operational observation window in which the deviation criterion is met before the waveform falls below the adopted amplitude threshold. The dependence of $t_{\mathrm{dev}}$ on the bump location is nonmonotonic. For example, with $A = 0.004$, $W = 1.25$, and the initial data, observation prescription, and thresholds adopted here, the condition $t_{\mathrm{peak}} \le t_{\mathrm{dev}} < t_{\mathrm{faint}}$ holds approximately for $10.6 < r_0/M < 21.3$. Changing the width can enlarge this interval or remove it over the sampled range. Thus, substantial QNM frequency shifts alone do not establish that the associated waveform differences emerge while the signal retains an appreciable relative amplitude. (Conversely, waveform deviations appearing before substantial signal decay do not by themselves imply sufficiently large frequency shifts to be resolved.)

The examples in § 1.3 demonstrate the value of assessing spectral shifts and waveform timing together. Potentially resolvable frequency shifts can coexist with waveform deviations emerging before the signal becomes too faint, and our matched-amplitude comparison shows that case (ii) can produce substantially larger shifts than case (i). This suggests that the bump effect may provide a route to amplify the weak quantum signatures expected in physically natural regimes of many LQBH models. Although case (ii) is phenomenological, it is admissible within the cLQBH framework and illustrates a generic response to localized potential deformations rather than an artifact of the particular profile. The bowl-shaped windows further suggest a selection mechanism that favors parameter choices balancing the magnitude of the departure from GR with its emergence time, rather than simply maximizing the frequency shift. The appearance of a similar window in the effective odd-parity calculation supports its relevance beyond the scalar examples without establishing universality. In turn, the cLQBH model provides a concrete BH setting for advancing bump physics from spectral characterization toward observationally oriented waveform analysis.

The observation window is a diagnostic of the calculated scalar response, and its boundaries depend on the stated deviation and faintness criteria. Likewise, our comparisons with reference accuracies motivated by LIGO observations and TianQin forecasts are exploratory and do not establish observational exclusions of the model. Quantitative gravitational-wave predictions require gravitational perturbation equations consistent with the cLQBH dynamics, extensions to rotating backgrounds, and analyses incorporating mode excitation, source properties, and detector noise. Alongside these developments, further investigations could extend the scalar analysis to nonminimal coupling, systematically search for exceptional points, and explore the possible connection between localized potential deformations and effective reflectivity [99]. The present results provide a framework for connecting localized holonomy corrections, QNM spectral changes, and the temporal development of waveform deviations, helping identify configurations that warrant more complete investigations.

Acknowledgments

C.Z. would like to thank Profs. WANG Peng, LIN Kai and YANG Runqiu, and Drs. GAN Wencong, MAI Zhanfeng and WU Liangbi for valuable discussions. A.W. and C.Z. also want to thank Prof. Shinji Mukohyama for suggestions and valuable discussions. C.Z. is supported in part by the National Natural Science Foundation of China (NSFC) under Grant No. 12205254, and the Science Foundation of China University of Petroleum, Beijing, under Grant No. 2462024BJRC005. A.W. is partially supported by the US National Science Foundation (NSF) through the grant: PHY-2308845.