In this study we investigate the quasinormal modes QNMs pseudospectrum and time evolution of a massive vector field around a quantum corrected black hole in deSitter spacetime. We start by parametrization and using orthonormal tetrads to get the effective potential. Methodologically we use the hyperboloidal framework together with discretizing the nonselfadjoint operator through ChebyshevGaussLabatto grid to attain the QNMs. We explore the parametric instability of QNMs caused by quantum correction cosmological constant and Proca mass and these three factors show very different influences on the QNMs migration flow. On the other hand we discuss the instability of QNMs with arbitraryshape perturbation and the effectiveness of numerical results through pseudospectrum. We use high frequency approximation to attain the expression of the time domain Green function and clarify the origin of two different stages in time evolution. Through numerical methods we confirm that no powerlaw late time tail is expected and the possible impact on time evolution caused by quantum correction is discussed.
['General relativity', 'alternative theories of gravity', 'General Relativity and Quantum Cosmology']
Quasinormal modes pseudospectrum and time evolution of Proca fields in a quantum OppenheimerSnyderde Sitter spacetime
{'The quasinormal modes, pseudospectrum and time evolution of Proca fields in quantum Oppenheimer-Snyder-de Sitter spacetime': "Shu Luo a ∗  \na School of Physics , University of Science and Technology of China, Hefei, Anhui 230026, China (Dated: September 13, 2024)  \nIn this study, we investigate the quasinormal modes, pseudospectrum and time evolution of a massive vector field around a quantum corrected black hole in de-Sitter spacetime. We start by parameterization and using orthonormal tetrads to get the effective potential. Methodologically we use the hyperboloidal framework together with discretizing the non-selfadjoint operator through Chebyshev-Gauss-Labatto grid to attain the QNMs. We explore the parametric instability of QNMs caused by quantum correction, cosmological constant and Proca mass, and these three factors show very different influences on the QNMs' migration flow. On the other hand, we discuss the instability of QNMs with arbitrary-shape perturbation and the effectiveness of numerical results through pseudospectrum. We use high frequency approximation to attain the expression of the time domain Green function and clarify the origin of two different stages in time evolution. Through numerical methods we confirm that no power-law late time tail is expected, and the possible impact on time evolution caused by quantum correction is discussed.", 'I. INTRODUCTION': "Recently, with the detection of gravitational waves(GW) from binary compact stars [1], together with the Event Horizon Telescope (EHT) revealing the optical characteristics of M87 center supermassive black hole [2], black holes, especially astronomical realistic black holes, can serve as a lab to test various theoretical predictions through Multi Messenger Astronomy, including quantum gravity and other modified gravity theories . These modifications mainly consists of two types: adding new elements to the right hand of the Einstein equation, i.e., supposing new source of energy-momentum tensor to tangle with gravitation, including dark matter and dark energy [3, 4]; modifying the left side of the equation by the new gravitational degrees of freedom in addition to the metric tensor. In many cases, the modification of GR can be described by a scalar-tensor theory of gravitation [5, 6]. Apart from the extensively-investigated scalar-tensor theory, there is another class of generalized vector-tensor theories, Einstein-Proca theories, based on the generalized Proca equation and higher order coupling of gravitation and Proca fields [7, 8], and there are investigations on relational spherically symmetric solutions [9], superradiant instabilities which offer restrict on the upper bounds of photons [10, 11], vector field quasinormal modes [12] and optical characteristics [13].  \nAmong these efforts, one important problem to be solved by new theories is spacetime singularities, characterized by infinite curvature or density, which has been a subject of great interest and curiosity in the fields of gravitation theory and relativistic astrophysics. The existence of singularities is unavoidable in generic gravitational collapses, which has been proven by Hawking and Penrose in the most general case. The presence of singularities poses profound challenges to our understanding of the universe within the context of classical general relativity, however, a probable candidate for quantum gravity theory, loop quantum gravity (LQG), provides a different view of this problem. It claims that through quantum geometry the singularity no matter in a black hole or the beginning of our universe can be solved [14, 15].  \nRecently, Lewandowski et al. have proposed a new quantum black hole model (qOS) within the framework of Loop Quantum Gravity (LQG) theory [16]. The exterior spacetime of this model is a suitably deformed Schwarzschild black hole. This quantum corrected black hole exhibits the same asymptotic behavior as the Schwarzschild black hole and it is stable against test scalar and vector fields by the analysis of QNMs [17]. Through the similar method, the solution describing a quantum black hole with a positive cosmological constant (qOS-dS) is obtained in [18].  \nAlthough having a lot of similarities to the common case, such a solution still has several observable properties that are different from those of the Schwarzschild-dS black hole, including its shadow and photon sphere, and by adding the cosmological constant, the authors in [18] focus on checking the SCC hypothesis in this spacetime. They put forward that SCC hypothesis will be destroyed as the black hole approaches the near-extremal limit. Additionally, in [18] some aspects about the quasinormal modes (QNMs) in the spacetime are also studied, but with few discussions on the influence of the perturbations of effective potential and the (in)stability of the QNMs, or quantitative description of the instability through pseudospectrum method.  \nThe studies of quasinormal modes (QNMs) mainly focus on characteristic behaviors of the modes under different types of spacetime. Starting from the 1970s, a lot of different spin fields' QNMs in a lot of kinds of D-dimensional spacetime have been deeply studied [19, 20]. Another center of investigation is the pseudospectrum and instability of the QNM, starting from the pioneering work of Nollert [21], which has shown that the QNM overtones are strongly unstable with their instabilities increasing with their damping [22, 23]. A very tiny parameter modification might sharply influence the behavior of QNMs, especially for the case of the comparison between a strictly massless field and very small but nonzero massive fields. It has been shown that  \nmassive fields own some modes that almost do not decay, corresponding to imaginary part very close to zero, which is called 'quasi-resonances' [24, 25]. This behavior very likely happens when the amplitude near the event horizon is even smaller than far from the black hole, which means there is no leak of energy from the system and there exist some modes of standing waves. Theoretical analyses also support purely real modes when there is a massive field [25].  \nWe mainly investigate the QNMs, their instability and evolution in the time domain of the quantum corrected black hole (qOS-dS) in this paper, and we choose the Proca field as the test field and explore their axial perturbation. This means we will combine the three important study center mentioned in the above paragraphs. And we will summarize the influence of these factors on QNMs by three parameters: the relative ratio of Cauchy horizon radius to outer horizon radius q , the ratio of outer horizon radius to cosmological horizon radius p and Proca mass m . The first reason of choosing this kind of field is to confirm the origin of the quasi-resonance phenomenon [24, 25] in a non asymptotic flat spacetime, and it will be shown that particularly for vector field the perturbation caused by changing Proca mass m can be analogized to continuous migration of the angular momentum number l , but m no longer plays the important role of deciding the asymptotic behavior, so no quasi-resonance is found. More practically, the detection on the Ringdown stage of the time domain (which has been known to closely dependant on QNMs) to ensure the upper limit of the quantum correction may be disturbed by the possible mass term of the photons, so it's important to clarify the influence caused by mass term.  \nThe first step of our study is the parameterization of the whole metric using the parameters p and q in the last paragraph, which can simplify the following investigations on the perturbations caused by quantum correction, excluding nonphysical range of parameters that do not accord to the premise discussed. As we need to construct a hyperboloidal coordinate to absorb the asymptotic boundary condition into non-divergent coordinate condition, compared to previous investigation [18] this parametric form seems to be more convenient. Then the governing equation of the axial perturbation of the Proca field in qOS-dS spacetime is obtained, together with using hyperboloidal framework to attain the QNMs. Moreover, in this process a possibility of losing the effectiveness of this method under the minimal gauge is discussed, which is of important physical meaning and offers a useful warning in later exploration.  \nWealso investigate the (in)stability of the spectrum from two aspects. On the one hand, we will use this solution as a correction to the Schwarzschild solution to discuss spectrum (in)stability of Schwarzschild black hole and the migration flow of the QNMs. In this way we will clarify the total tendency of the frequency domain when some parameters migrate. And some important characteristics under quantum correction is revealed. On the other hand, we will use the pseudospectrum method to study the relative stability of the QNMs of different overtones under arbitrary perturbation of a fixed amplitude, which has never been done before for this solution.  \nAfter all these discussions are finished, the evolution in the time domain is also discussed in detail. We mainly focus on the different stages of the evolution. The lasting time of the first stage, Precursor, is influenced by parameters. The second stage, Ringdown, also show some important signals affected by quantum correction. About the third stage, as mentioned in [26], the late-time tails can reflect some essential properties of black holes, such as the no-hair theorem and the instability of Cauchy horizons, and here under a non-zero cosmological constant it strictly obeys exponentially decaying rule. Also, a detail about the distribution of the energy on the spatial domain is incidentally mentioned afterwards.  \nThe paper is organized as follows. In the next section we finish the parameterization of the qOS-dS spacetime. Then in Sec.III, we present the equation of motion of Proca field in qOS-dS spacetime. In Sec.IV we construct the hyperboloidal coordinate to calculate QNMs numerically. In Sec.V we discuss the instability of the QNMs, first by parametric migration in Sec.V A and then by pseudospectrum in Sec.V B. We then attain the high-frequency approximation of the Green function and numerically calculate the evolution in time domain in Sec.VI. Sec.VII is the conclusions and discussion.", 'II. PARAMETERIZATION OF THE QOS-DS SOLUTION': 'We first briefly discuss the physical image of a qOS-dS model. Through the GR, it has been revealed that a collapsing star can be described by the Oppenheimer-Snyder model, which predicts that a singularity in the final evolutionary stage of the star is unavoidable. However, the loop quantum theory (LQG) predicts a different picture of the collapsing and bouncing matters in a certain dust ball of a time-changing radius never smaller than a certain r b . According to the Ashtekar-Pawlowski-Singh (APS) model, the inner region of a such dust ball can be described by an non-globally static metric in the same form of the RobertsonWalker metric with k = 0 , while the exterior spacetime is a pseudo-static one, which has a form close to Schiwarzschild spacetime. The only difference lies in the dynamical equation ruling the evolution of the scale parameter, or, equivalently, the radius of the dust ball.  \nAnd the mere difference between the qOS-dS and qOS model is the non-vanishing cosmological constant, which contributes to the dynamical evolution of the scale factor a ( τ ) , and the expression of the exterior spacetime is also modified. Specifically speaking, the metric has the form of  \nd s 2 in = -d τ 2 + a ( τ ) 2 [d˜ r 2 + ˜ r 2 (d θ 2 +sin 2 θ d ϕ 2 )] , (2.1)  \ninside the dust bull surface, in which the scale factor is governed by (we have set gravity constant G = 1 )  \nH = ( ˙ a a ) 2 = 8 π 3 ρ (1 -ρ ρ c ) + Λ 3 , ρ = M 4 π 3 a 3 ˜ r 3 0 , (2.2)  \nin which ρ c = √ 3 / 32 π 2 γ 3 l 2 p , l p = √ ℏ is the Plank length, γ being the Immirzi parameter, Λ being the cosmological constant and M standing for the mass of the dust ball. The coordinate ˜ r only makes sense in the area 0 < ˜ r < ˜ r 0 .  \nFigure 1: Penrose diagram for the maximum extension of the spacetime with three horizons in the exterior area of the dust bull. However, only the areas A ext , B ext and C ext and D are real and the metric (2.3) can make sense.  \n<!-- image -->  \nBut in our paper, we focus on the exterior spacetime of the qOS-dS black hole whose metric is given by [18]:  \nd s 2 ex = -f ( r )d t 2 + d r 2 f ( r ) + r 2 (d θ 2 +sin 2 θ d ϕ 2 ) , (2.3)  \nwhere the metric function f ( r ) reads  \nf ( r ) = 1 -2 M r -Λ 3 r 2 + αM 2 r 4 ( 1 + Λ r 3 6 M ) 2 , (2.4)  \nwith the positive parameter α = 16 √ 3 πγ 3 l 2 p .  \nThe most realistic case is when there is three horizons, corresponding to three positive roots of the metric. If we set r -, r + and r 0 in response to the inner, outer and cosmological horizon respectively, then one special case is that when Λ = 0 and the spacetime is asymptotically flat, then there are only two horizons expected, under which following the method in [27] and setting 0 < q 2 = r -/r + < 1 one can easily get the relation  \nα M 2 = 16 q 6 (1 + q 2 + q 4 ) 3 (1 + q 2 + q 4 + q 6 ) 4 , (2.5)  \nand  \nr + = 2 M 1 + q 2 + q 4 1 + q 2 + q 4 + q 6 , r -= 2 Mq 2 1 + q 2 + q 4 1 + q 2 + q 4 + q 6 , (2.6)  \nand the more general case is that there is also a cosmological horizon, for which through similar strategy, that is, setting 0 < q 2 = r -/r + < 1 and 0 < p 2 = r + /r 0 < 1 , then the metric can be rewritten as  \nf ( r ) = ( α Λ 2 36 -Λ 3 ) r 2 ( 1 -r 0 r )( 1 -p 2 r 0 r )( 1 -p 2 q 2 r 0 r )( 1 + dr 0 r + br 2 0 r 2 + cr 3 0 r 3 ) , (2.7)  \nwhere  \nd = 1 + p 2 + q 2 p 2 , b = p 2 q 2 (1 + p 2 + q 2 p 2 )(1 + q 2 + q 2 p 2 ) 1 + q 2 + q 4 + p 2 q 2 + p 2 q 4 + p 4 q 4 , c = p 4 q 4 (1 + p 2 + q 2 p 2 ) 1 + q 2 + q 4 + p 2 q 2 + p 2 q 4 + p 4 q 4 , (2.8)  \nin this way one is easy to find the relation of parameters including  \nαM 2 r -4 0 = ˜ a ˜ c , ( 2 -α Λ 3 ) M r 0 = ˜ b ˜ c , ( Λ 3 -α Λ 2 36 ) r 2 0 = ˜ c , (2.9)  \nwhere  \n˜ a = p 8 q 6 (1 + p 2 + q 2 p 2 ) 1 + q 2 + q 4 + p 2 q 2 + p 2 q 4 + p 4 q 4 , (2.10)  \n˜ b = p 2 (1 + p 2 )(1 + q 2 ) 1 + q 4 + p 4 q 4 + p 6 q 6 + p 2 ( q 2 + q 6 ) 1 + q 2 + q 4 + p 2 q 2 + p 2 q 4 + p 4 q 4 , (2.11)  \n˜ c = 1 + q 2 + q 4 + p 2 q 2 + p 2 q 4 + p 4 q 4 (1 + p 2 + p 4 )(1 + q 2 + q 4 )(1 + p 2 q 2 + p 4 q 4 ) . (2.12)  \nBy combing the relations in Eqs.(2.9), one is able to get more explicit relations between given parameters:  \nΛ α = 6 ( 1 -√ 1 1 + 4˜ a/ ˜ b 2 ) , α M 2 = 16˜ a ˜ c 3 ( ˜ b 2 +4˜ a ) 2 , r 0 M = 2 ˜ c √ 1 ˜ b 2 +4˜ a . (2.13)  \nFrom Fig.2, it is shown that the largest value of Λ α is 6 -4 √ 2 , if and only if p 2 = q 2 = 1 , i.e., the three horizons merge. At the same time, we have α/M 2 = 9 / 4 and r 0 /M = 3 √ 2 / 2 as p 2 = q 2 = 1 . Through comparing the top right panel and the bottom right panel, one is able to find as long as the black hole mass M is fixed, p mostly relies on parameter Λ , while q mostly relies on the parameter α . That is a natural conclusion considering the physical meaning of these two factors. The main advantage of this parameterization is the convenience in the following numerical calculation, in which extreme case is easier to study.  \nFor the details of this spacetime, we can effectively adopt the semi-classical method and suppose the metric of the spacetime still obey the Einstein equation with nothing but a nonzero quantum energy-momentum tensor [16]. Now one can easily calculate the non-zero components of the energy momentum tensor T µν as  \nT 00 = (36 αM 2 +12Λ r 6 -α Λ 2 r 6 )[ α (6 M +Λ r 3 ) 2 -12 r 3 (6 M -3 r +Λ r 3 )] 432 r 10 +Λ ( 1 -2 M r -Λ 3 r 2 + αM 2 r 4 ( 1 + Λ r 3 6 M ) 2 ) , (2.14)  \nT 11 = 3Λ( α Λ -12) r 6 -108 αM 2 α (6 Mr + αr 4 ) 2 -12 r 5 (6 M -3 r +Λ r 3 ) + Λ 1 -2 M r -Λ 3 r 2 + αM 2 r 4 ( 1 + Λ r 3 6 M ) 2 , (2.15)  \nT 22 = 6 αM 2 r 4 + Λ 2 12 α, T 33 = [72 αM 2 +Λ 2 αr 6 ]sin 2 θ 12 r 4 +Λ r 2 sin 2 θ , (2.16)  \nand some of important scalars include  \nR = 4Λ -α Λ 2 3 -6 M 2 α r 6 , (2.17)  \nR ab R ab = 90 M 4 α 2 r 12 + M 2 α Λ( α Λ -12) r 6 + Λ 2 ( α Λ -12) 2 36 , (2.18)  \nR abcd R abcd = 468 M 4 α 2 r 12 + 40 M 3 α ( α Λ -6) r 9 + 2 M 2 ( α 2 Λ 2 -12 α Λ+24) r 6 + Λ 2 ( α Λ -12) 2 54 , (2.19)  \nFigure 2: The relation between some given parameters and the relative ratio of the three horizon radius.  \n<!-- image -->', 'III. AXIAL PROCA FIELD PERTURBATION OF QOS-DS MODEL': "In this section, we will study the perturbation of the Proca field in the qOS-dS spacetime. The premise of our study is the semi-classical approximation, that is, to regard the Proca field as a classical field on a spacetime metric influenced by quantum effect. Following the method proposed in [28], for the spacetime described by  \nd s 2 = -f ( r )d t 2 + d r 2 g ( r ) + r 2 (d θ 2 +sin 2 θ d ϕ 2 ) , (3.1)  \nwe study the linear Einstein-Proca equation in an orthonormal frames (in the following calculation we use f to stand for f ( r ) and g to stand for g ( r ) ):  \ne (0) µ = ( -√ f, 0 , 0 , 0) , e (1) µ = (0 , 1 √ g , 0 , 0) , e (2) µ = (0 , 0 , r, 0) , e (3) µ = (0 , 0 , 0 , r sin θ ) . (3.2)  \nor, equivalently,  \ne µ (0) = ( 1 √ f , 0 , 0 , 0) , e µ (1) = (0 , √ g, 0 , 0) , e µ (2) = (0 , 0 , 1 r , 0) , e µ (3) = (0 , 0 , 0 , 1 r sin θ ) . (3.3)  \nin which later calculation shows obvious simplification in the equation's expression, as there is e µ ( a ) e ( b ) µ = η ( a )( b ) . Now we turn to the study of the axial perturbation of Proca field. 'Axial' means the rotation of the field around a certain axis, which contributes to an only nonzero component of vector potential in A 3 , and this component can contribute to nonzero F 03 , F 13 and F 23 . If the original A 0 is nonzero, then this perturbation can also contribute to the non diagonal metric components including g 03 , g 13 , g 23 but no other components. It can be proved that this kind of perturbation is independent of another perturbation called polar perturbation [28]. In the frame introduced earlier, the Proca equation has the form of  \nη ( m )( n ) F ( a )( m ) | ( n ) + m 2 A ( a ) = 0 , F [( a )( m ) | ( n )] = 0 , (3.4)  \nHowever, as the specific expression of this equation requires the calculation of the rotation coefficients, we firstly adopt the common Proca equation in the coordinates:  \n∇ µ F µν + m 2 A ν = 0 , ∇ [ µ F νρ ] = 0 , (3.5)  \nor, explicitly, for axial perturbation:  \n√ g f 1 r 2 ∂ ∂r ( √ f g r 2 sin θF 13 ) + ∂ ∂θ (sin θF 23 ) + ∂ ∂t (sin θF 03 ) + m 2 A 3 = 0 , (3.6)  \ntogether with  \nConsidering the relation  \nthus all we need are  \nm 2 √ f g rA (3) + r √ g F (0)(3) ,t +( r √ fF (3)(1) ) ,r + √ f g F (3)(2) ,θ = 0 , (3.9)  \n( r √ fF (0)(3) sin θ ) ,r + r √ g F (3)(1) ,t sin θ = 0 , r √ f ( F (0)(3) sin θ ) ,θ + r 2 F (3)(2) ,t sin θ = 0 , (3.10)  \ninsert (3.10) ( ?? ) into (3.9) and note that A (3) ,t = √ fF (0)(3) we get  \nr [ Ω 2 √ g -m 2 f √ g ] F (0)(3) + f √ gr ∂ ∂θ [sin θ ∂ ∂θ ( F (0)(3) sin θ )] + ∂ ∂r [ √ fg ∂ ∂r ( r √ fF (0)(3) )] = 0 . (3.11)  \nhere we have supposed a feature frequency ω to replace the derivative to t . Then by variable separating we get  \nd d r [ √ fg d d r ( r √ fB )] -l ( l +1) fB √ gr + r ( ω 2 √ g -m 2 f √ g ) B = 0 . (3.12)  \nwhere B is the radial part of the whole perturbation, l is the angular momentum number.", 'IV. NUMERICAL APPROACH TO CALCULATE QNMS OF AXIAL PROCA FIELD PERTURBATIONS IN QOS-DS BLACK HOLE': "Firstly, we should make it clear that for a non vanishing Λ the region we study is between the outer horizon ( r + ) and cosmological horizon ( r 0 ) of the qOS model. Under this basic assumptions, there is always V ( r ) → 0 no matter for r → r + or r → r 0 , so both sides are null boundary condition, and more specifically, there is Φ ∼ e iω ( t -r ∗ ) near r = r 0 and Φ ∼ e iω ( t + r ∗ ) near r = r + , which accords to conventional definition for QNMs, even though here the field is massive. This offers the central prerequisite for using hyperboloidal approach, which are fundamentally revealed in [ ? ]. The hyperboloidal approach is able to absorb the boundary asymptotic condition into the non-divergence condition of the differential operator near the boundary, which is naturally satisfied in discretizing and matrix approximation. However, for Λ = 0 and an asymptotically flat spacetime, the boundary condition is time-like for massive fields, which cause the traditional hyperbolodial methods in calculating QNMs lose effectiveness. Here we avoid to discuss this condition. For the special case f = g , we set  \nΦ = -r √ fB 2 √ n , d r ∗ = d r f , (4.1)  \nand insert into (3.12), we get the wave equation  \n[ ∂ 2 ∂r 2 ∗ + ω 2 -f r ( l ( l +1) r + m 2 r )] Φ = 0 , (4.2)  \n∂ r F 03 + ∂ t F 31 = 0 , ∂ θ F 03 + ∂ t F 32 = 0 , (3.7)  \nF µν = F ( a )( b ) e µ ( a ) e ν ( b ) = η ( c )( a ) η ( d )( b ) F ( c )( d ) e µ ( a ) e ν ( b ) , (3.8)  \nWhen m = 0 , the potential is exactly the Regge-Wheeler potential for s = 1 . As the space region we study is r + < r < r 0 , using the minimal gauge of hyperboloidal we transform the coordinate to  \nr ∗ = 2 r + h 1 ( σ ) , t = 2 r + ( τ -h 2 ( σ )) , (4.3)  \nusing the minimal gauge, the explicit form of the transformation can be written as  \nσ = r + r 0 ( r 0 -r + ) r -r + r 0 -r + , d τ = d t -d r ∗ 2 r + + βr + d r r ( r -r + ) , (4.4)  \nwhere  \nβ = -1 -p 2 f ' (1) , (4.5)  \nin which p 2 = r + /r 0 and f is written as a function of σ . Fig.3 shows the effect of this transformation: the original infinity in t -r domain is drawn closer and divided into several τ = constant hypersurfaces.  \nFigure 3: The τ = constant hypersurfaces in t -σ diagram under the parameters p = q = 0 . 5 . It's shown that while τ remains constant t approaches infinity on both sides.  \n<!-- image -->  \nFollowing the method in [27, 29], we set Ψ = ∂ τ Φ , then in ( σ, τ ) coordinate there is  \niω [ Φ Ψ ] = ∂ τ [ Φ Ψ ] = [ O 1 L 1 L 2 ][ Φ Ψ ] ≡ L [ Φ Ψ ] , (4.6)  \nwhere  \nand  \nL 1 = 1 w 1 ( σ ) ( ∂ σ ( ν ( σ ) ∂ σ ) -Q ( σ )) , (4.7)  \nL 2 = 1 w ( σ ) (2 γ 1 ( σ ) ∂ σ + ∂ σ γ ( σ )) , (4.8)  \nν ( σ ) = 2[(1 -p 2 ) σ + p 2 ] 2 f 1 -p 2 , Q ( σ ) = 2 r 2 + (1 -p 2 ) V [(1 -p 2 ) σ + p 2 ] 2 f , (4.9)  \nw ( σ ) = 2 β 1 -σ -2 β 2 [(1 -p 2 ) σ + p 2 ] 2 f (1 -σ ) 2 (1 -p 2 ) , γ ( σ ) = 1 -2 βf [(1 -p 2 ) σ + p 2 ] 2 (1 -p 2 )(1 -σ ) . (4.10)  \nand here ω actually stands for dimensionless 2 r + ω .  \nHowever, the effectiveness of the method in Eqs.(4.9) and Eqs.(4.10) relies on the non-divergence of the elements in the calculating matrix in Eq.(4.6). Although the setting of the coefficient in Eq.(4.5) have already precisely guaranteed that w ( σ ) , γ ( σ ) and their reciprocals are regular on the two edges and the vanishing of the effective potential on the edges guaranteed the non-  \nergence of Q ( σ ) , there is still a residual danger that under certain parameters whether w ( σ ) could have a root between σ = 0 and σ = 1 . Unfortunately for qOS-dS model there really exists such a forbidden region of parameters p and q where our hyperboloidal coordinates construction loses effectiveness, which is shown in Fig.4.  \nFigure 4: The p -q diagram for the model, where the orange region corresponds to the permitted region for our method while the blue region corresponds to the forbidden region.  \n<!-- image -->  \nSo what happens when the parameters are in the forbidden region? Actually one is easy to confirm that according to the general transformation rule in Eqs.(4.3) the hypersurface τ = constant is marked by its normal vector (d τ ) a , and there is  \ng ab (d τ ) a (d τ ) b = -1 (2 r + ) 2 f h ' 2 1 -h ' 2 2 h ' 2 1 -h ' 2 2 + h ' 2 1 h ' 2 2 , (4.11)  \nin which h ' 1 and h ' 2 stands for the derivatives to coordinate σ . Meanwhile there is  \nw ( σ ) = h ' 2 1 -h ' 2 2 | h ' 1 | , (4.12)  \nUsually there is g ab (d τ ) a (d τ ) b < 0 , which means the τ = constant surface is a space-like one and accords to the demand, however, this is not always true. It might be guaranteed on one side (in our method σ = 0 ), when being violated on the other, which, according to the continuity, must lead to one root of g ττ as well as the function w ( σ ) between σ = 0 and σ = 1 . Under this condition our coordinate hypersurface turns out to be space-like only on partial σ domain but null and even time-like on the others. And in the previous studies although the asymptotic behavior of τ coordinate line on the two horizons is emphasized for many times, this fact is often neglected. Although this is a meaningful warning for coordinate constructing, here for convenience we just avoid talking about the forbidden region or using more complex hyperboloidal gauge.  \nTo calculate the QNMs, we firstly transfer the variable range from (0 , 1) to ( -1 , 1) :  \nx = 2 σ -1 , (4.13)  \nwe sample N +1 points from a Chebyshev-Gauss-Labatto grid in which the points are the extrema of Chebyshev polynomial on the interval ( -1 , 1) , i.e.,  \nx j = cos ( jπ N ) , j = 0 , 1 , · · · , N , (4.14)  \nand follow [29] we can discretize the original derivative matrix equation (4.6) to  \niω [ Φ Ψ ] = L [ Φ Ψ ] = [ O I N L (2) 1 D 2 N + L (1) 1 D N + V L (1) 2 D N + L (0) 2 ][ Φ Ψ ] , (4.15)  \nwhere  \nΦ = (Φ( x 0 ) , Φ( x 1 ) , · · · , Φ( x N -1 ) , Φ( x N )) T , (4.16)  \nand  \nand  \nΨ = (Ψ( x 0 ) , Ψ( x 1 ) , · · · , Ψ( x N -1 ) , Ψ( x N )) T , (4.17)  \n( D N ) 00 = 2 N 2 +1 6 , ( D N ) NN = -2 N 2 +1 6 ,  \n( D N ) jj = -x j 2(1 -x 2 j ) , j = 1 · · · , N -1 ,  \n̸  \n( D N ) ij = c i c j ( -1) i + j ( x i -x j ) , i = j , i, j = 0 , · · · , N , (4.18)  \n(4.19)  \nc i = { 2 , i = 0 or N , 1 , otherwise . (4.20)  \n( L (2) 1 ) ij = 4 α ( x i +1 2 ) δ ij w ( x i +1 2 ) , ( L (1) 1 ) ij = 4 ∂ x α ( x i +1 2 ) δ ij w ( x i +1 2 ) , (4.21)  \n( L (1) 2 ) ij = 4 γ ( x i +1 2 ) δ ij w ( x i +1 2 ) , ( L (0) 2 ) ij = 2 ∂ x γ ( x i +1 2 ) δ ij w ( x i +1 2 ) , (4.22)  \n( V ) ij = Q ( x i +1 2 ) δ ij w ( x i +1 2 ) . (4.23)", 'A. Parameter Perturbation': "About the stability under a certain perturbation, we introduce the migration susceptibility defined as  \nχ n ≡ lim ϵ → 0 δω n ϵ , (5.1)  \nto quantify the relative instability of a certain QNM, where n stands for the overtone number and δω n refers to the change between unperturbed QNMs and the perturbed ones. It should be noted that the sort of overtone is still based on the absolute value of the imaginary part. In this definition, the fundamental mode is considered as the mode with the smallest imaginary part. For the most general discussion, we replace r ∗ by x , and suppose the original potential to be V ( x ) with the migration to be ϵV b ( x ) . The definition of Jost solutions is those solutions of Eq.(4.2) satisfying the conditions [30]  \nϕ ω ( x ) → e iωx [ 1 + O ( 1 x ) ] , x →-∞ (5.2)  \nψ ω ( x ) → e -iωx [ 1 + O ( 1 x ) ] , x → + ∞ (5.3)  \nand as there are only two linear independent solutions to wave equation, there is the relation  \nϕ ω ( x ) = a ( ω ) ψ -ω ( x ) + b ( ω ) ψ ω ( x ) , (5.4)  \nand through this definition it's easy to see all the QNMs are the roots of the a ( ω ) , so considering the perturbation from a (0) ( ω ) to a ( ϵ ) ( ω ) caused by ϵV b ( x ) , there is the relation  \n0 = a (0) ( ω n ) = a ( ϵ ) ( ω ( ϵ ) n ) ≈ a ( ϵ ) ( ω n ) + a ( ϵ ) ' ( ω n ) δω n , (5.5)  \nso there is  \nχ n = -lim ϵ → 0 [ 1 ϵ a ( ϵ ) ( ω n ) a (0) ' ( ω n ) ] , (5.6)  \nand from the definition in (5.2) (5.3), one is able to solve the explicit form of the coefficient a ( ω ) :  \na ( ω ) = 1 -∫ ∞ -∞ e iωy 2 iω V ( y ) ϕ ω ( y )d y , (5.7)  \nand the numerator of χ n can be approximated as  \na ( ϵ ) ( ω n ) ≈ -ϵ ∫ ∞ -∞ e iω n y 2 iω n V b ( y ) ϕ ω n ( y )d y , (5.8)  \nand is usually a regular number compared to the order ϵ as long as V b is also a short-ranged and positive potential so it does not contribute much to either stability or the instability. About the denominator, however, things are quite different, as there is  \na (0) ' ( ω n ) = -∫ ∞ -∞ e iω n y 2 iω n V ( y ) [ iω n y -1 ω n ϕ ω n ( y ) + δϕ ω n ( y ) δω n ]d y , (5.9)  \nand if this denominator turns to be 0, then there is a so-called Type I instability [30], in which χ n under a perturbation turns to be singular, and it's clear that such singularity has nothing to do with the specific form of perturbation.  \nNow that we have constructed the hyperboloidal framework in Sec.IV, we can use the calculation results of the first overtones to study the instability and their response to parameter migration. We take N = 48 as the dimension of the discrete matrix L and machine precision to be 120 and we confirm that at least for the first four overtones the results are convergent as N increases.  \nFrom the top right and bottom right panel of Fig.2 one can confirm that when black hole mass M is fixed, the cosmological constant Λ mainly relies on p , while α mainly relies on q . For specific perturbations we firstly study the corrections caused by quantum effect itself to the original Schiwarzschild-dS QNMs, under different overtone number n and angular momentum number l . So we change parameter q under a constant parameter p , and when q reduces to 0 it's just the case for SchiwarzschilddS spacetime. Also, one can study the perturbation caused by cosmological constant: when q is fixed and p changes from 0 to 1, and the results are shown in Fig.5.  \nThe general rules are as follows: the QNMs of the Proca field are more unstable under parametric perturbations caused by p than by q . It's rather easy to understand the reason: the change in p finally ends up in a quite different result in the asymptotic behavior of the Proca field: when p =0 the Proca behaves as the really massive field as in the flat spacetime. For p or q , when they approach 1 together, the three horizons merge and the effective potential diminishes, so the QNMs reduces to near 0 and is reflected by Fig.5. In the migrating process there are several points for these first overtones to behave Type I instability where there are obvious discontinuity on some of these curves, and this is more likely to happen for higher overtones. By comparing the shapes of the migration curves in Fig.5 with some previous works of the pseudospectrum of QNMs [27], it's very clear that even for the most stable fundamental modes, they still migrate in the unrestricted lines. And this originates from the fact that the operator represented by the discontinuous matrix L is non self-adjoint, and there is no guarantee of the stability of its eigenvalues under perturbation.  \nFigure 5: The migration of the first 3 or 4 overtones under the parameter m,l fixed and simply changing p from 0 to 1 (left) or simply changing q from 0 to 0.8 (right). Here and in the following calculation q does not range from 0 to 1 to avoid being captured in the forbidden region.  \n<!-- image -->  \nAnother aspect to study is the mass constant m , which, according to the expression of the effective potential, can simply  \nreduce the possible difference caused by different l up to a common ratio for all QNMs, so we don't need to change l any longer. In Fig.6, the QNMs migration of first overtones under the change in m is shown, and it is clear that m mostly contributes to a increase in the real part of the QNMs but does nothing to the imaginary part, so the increase of the oscillating rate is expected under the growth in the Proca mass, but no change in the decay rate. The rule for massive fields, the first overtone of the QNMs becoming very small in its imaginary part as m increases, which, according to a lot of previous researches, are called quasiresonance, is no longer observed. We find m to have (at least for real part) influence nearly opposite from p and q , which warns that there may be a possible 'mass concealment' effect.  \nAfter attaining the general migration of the QNMs, it's meaningful to study the specific migration of the real parts and imaginary part of QNMs separately. Some results are shown in Fig.7 and Fig.8. In the first figure of Fig.8, a phenomenon of overtone replacement is observed, as when p approaches 0, Type I instability occurs and continuous curves of overtone migration are disrupted by suddenly occurring new branches. This phenomenon rapidly lowers the total imaginary part of QNMs, with greater possibility to move to quasi-resonance, which indicates the fundamental QNM ending up in imaginary part very close to 0. This, together with the result of changing m but keeping p a constant, once again illustrates the asymptotic behavior of the field is more important than the mass itself.  \nFigure 6: The migration of the first 3 or 4 overtones under the parameter l, p, q fixed and simply changing m from 0 to 1.  \n<!-- image -->  \nFigure 7: The change of the real part of the first 3 overtones under the parameter l, m fixed and simply changing p from 0 to 1 (left) or simply changing q from 0 to 0.8(right).  \n<!-- image -->", 'B. Pseudospectrum': "Given ϵ > 0 , the ϵ -pseudospectrum σ ϵ ( A ) of an operator A is defined as [27]  \nσ ϵ ( A ) = { z ∈ C : ∥ R A ( z ) ∥ = ∥ ( z I -A ) -1 ∥ > 1 /ϵ } , (5.10)  \nwhere R A ( z ) is called the resolvent operator. This definition is the most suitable for visualizing the pseudospectrum. In the limit ϵ → 0 , the set σ ϵ ( A ) reduces to the spectrum set σ ( A ) , whose elements are the spectrum ω n , in which n is the overtone number. The quantity ϵ serves as a measure of the 'proximity' between points in σ ϵ ( A ) and the spectrum ω n , offering a clear interpretation of perturbations to the underlying operator. Therefore, the shape and size of the ϵ -pseudospectrum regions quantify  \nFigure 8: The change of the imaginary part of the first 4 overtones under the parameter l, m fixed and simply changing p from 0 to 1 (left) or simply changing q from 0 to 0.8(right).  \n<!-- image -->  \nthe spectrum (in)stability of the operator A . If the contour lines of the pseudospectrum forms a concentric circle of the spectrum, then the operator considered is spectrally stable. On the contrary, if this structure does not occur, then this operator is spectrally unstable. It is known that with regard to the QNMs, the spectrum is generally unstable since the black hole is actually a leaky system [31, 32].  \nAbout addressing the construction of pseudospectrum in a numerical approach, we can use Chebyshev differentiation matrices and Chebyshev-Lobatto N -point grids to produce L matrix approximates. Once the operator is discretized, the construction of the pseudospectrum requires the evaluation of matrix norms. A standard practical choice involves the matrix norm induced from the 2-norm in the vector space, which leads to the following rewriting [29]:  \nσ ϵ 2 ( L ) = { z ∈ C : σ min ( z I -L ) < ϵ } , (5.11)  \nwhere σ min ( M ) stands for the smallest singular value of the matrix M . The most significant advantage to choose the 2-norm is the convenience and exemption from constructing 'energy norm', which is easy to loose effectiveness for systems with more than 1 degree of freedom, or the energy is not positive definite, such as in Kerr spacetime.  \nIn Fig.9 and Fig.10, the contour map of the pseudospectrum is shown, and each contour line corresponds to a constant σ min ( z I -L ) on the z plane. As the imaginary part increases, the instability of QNMs get more and more severe. The nonHermitian of the operator L is reflected by those unclosed contour lines, although from the right panel of Fig.9 one can still see very tiny enclosed circles around a QNM. These circles offer a important index of whether a numerical result of a QNM is believable: if the calculation error is found to be even larger than the outermost circle's ϵ value, than the result would be non-convergent. Also, the contour lines show pretty good accordance to the unclosed migration flow in Sec.V A, but here we can even extend the meaning of a perturbation to a arbitrary shape, given a certain amplitude. We can therefore conclude that the results shown in Fig.5 is not coincident.  \nNow that we have got the characteristics of the perturbation of QNMs in the frequency domain, we naturally want to get more knowledge about the evolution in the time domain as well, which will be illustrated in the next section.", 'VI. EVOLUTION IN THE TIME DOMAIN': "We first need to get the expression of the Green function in the frequency domain. By introducing s = iω , and remaining the definition of coefficients for a specific ω in Eq.(5.4) and Eq.(5.2), we have  \nG ( s, x, x ' ) = 1 2 sa ( s ) ( θ ( x -x ' ) ϕ s ( x ' ) ψ s ( x ) + θ ( x ' -x ) ϕ s ( x ) ψ s ( x ' )) , (6.1)  \nand then in time domain we can use the inverse Laplace transformation:  \nG ( t, x ' , x ) = 1 2 πi ∫ s 0 + i ∞ s 0 -i ∞ d s e st G ( s, x, x ' ) , (6.2)  \nwhere s 0 is a positive number. Through residue theorem we can transfer this to three parts, so in general, as illustrated in [30], there should be three main stages of the evolution: 1. Precursor, mainly coming from integration along the large semicircle on the half plane with | s | → ∞ . 2. Ringdown, mainly coming from the poles of a ( s ) , which is related to QNMs, (actually related to greybody factor which is stable under perturbation [33]). 3. Power-law tail, mainly coming from branch cut on the negative  \nFigure 9: The ϵ -pseudospectrum of the quantum corrected black hole with q = p = 0 . 5 , m = 10 -3 and l = 1 . The solid contour lines correspond to various choices of log 10 ϵ . The left panel represents the zoomed-out view of the ϵ -pseudospectrum. The right panel represents the zoomed-in view of the ϵ -pseudospectrum around the fundamental mode n = 0 .  \n<!-- image -->  \nFigure 10: The ϵ -pseudospectrum of the quantum corrected black hole with q = 0 . 1 , p = 0 . 8 , m = 10 -3 and l = 1 . The solid contour lines correspond to various choices of log 10 ϵ . The left panel represents the zoomed-out view of the ϵ -pseudospectrum. The right panel represents the zoomed-in view of the ϵ -pseudospectrum around the fundamental mode n = 0 .  \n<!-- image -->", 'real axis (NRA).': "However, as in the qOS-dS spacetime, the effective potential decays exponentially, there is actually no branch cut on the NRA. So there is no third stage at all, which means no power-law tail is expected. For convenience we focus on the case when the initial data mainly lies on positions far in the future of the peak of the potential. Actually later numerical calculation shows only little difference for the different observer's position in the main part of time, so we suppose the observer being in infinite future as well, but on the left side of the source. Now we only need to consider the asymptotic behavior of the Green function G ( t, x ' , x ) as x ' →∞ , x →∞ and x < x ' , where x ' is the position of the source and x is the position of the observer. So now through the definition in Eq.(5.4) we have  \nG ( s, x ' , x ) = 1 2 s [ ψ -s ( x ) ψ s ( x ' ) + b ( s ) a ( s ) ψ s ( x ) ψ s ( x ' ) ] , x, x ' →∞ , (6.3)  \nso we naturally divide Eq.(6.2) into two parts:  \nG (1) ( t, x ' , x ) = L -1 [ 1 2 s ψ -s ( x ) ψ s ( x ' ) ] , G (2) ( t, x ' , x ) = L -1 [ b ( s ) 2 sa ( s ) ψ s ( x ) ψ s ( x ' ) ] , (6.4)  \nFor high frequency limit | s | → ∞ one can suppose [34]  \nψ s ( x ) = e -sx [ 1 -g 1 ( x ) s + g 2 ( x ) s 2 ] + O ( 1 s 3 ) , x →∞ , (6.5)  \nand by inserting this into the wave equation we can solve g 1 ( x ) and g 2 ( x ) order-by-order in s :  \ng 1 ( x ( r )) = 1 2 [ l ( l +1) rr 0 + m 2 ] ( r -r 0 ) , (6.6)  \ng 2 ( x ( r )) = 1 8 ( r -r 0 ) 2 [ l ( l +1) r 0 r + m 2 ] 2 + 1 4 [ l ( l +1) r 2 + m 2 ] f ( r ) , (6.7)  \nin which we have guarantee that as r → r 0 both g 1 ( x ) and g 2 ( x ) reduce to 0. Simply replacing s by -s and r 0 with r + we get the expression for ϕ s ( x ) as x →-∞ . To get the expression for a ( s ) we need to approximate when x ≈ 0  \n2 sa ( s ) ≈ ψ s (0) d ϕ s (0) d x -ϕ s (0) d ψ s (0) d x = 2 s +( r 0 -r + ) [ l ( l +1) r + r 0 + m 2 ] + O ( 1 s ) , (6.8)  \n2 sb ( s ) ≈ ϕ s (0) d ψ -s (0) d x -ψ -s (0) d ϕ s (0) d x = r 0 -r + 2 s [ l ( l +1) r + r 0 + m 2 ] + O ( 1 s 2 ) , (6.9)  \nG (1) ( s, x ' , x ) ≈ e s ( x -x ' ) 2 s , G (2) ( s, x ' , x ) ≈ ( r 0 -r + ) 8 s 3 [ l ( l +1) r + r 0 + m 2 ] e -s ( x + x ' ) , (6.10)  \nand now we divide the whole time domain into three parts. First of all, when t < x ' -x , then for both G (1) ( t, x ' , x ) and G (2) ( t, x ' , x ) , we need to choose the large semicircle on the right-half plane so both semicircle contour integrals are 0 according to Jordan Theorem, with no residue in their contours as well. Now G ( t, x ' , x ) = 0 . This result is natural considering the not-so-good geometric optical approximation made in Eq.(6.5) when x ' is still finite.  \nWhen x ' -x < t < x ' + x , G (2) ( t, x ' , x ) is still 0 while for G (1) ( t, x ' , x ) we have to choose the semicircle on the left-half plane, and now the total contribution is  \nG (1) ( t, x ' , x ) = Res [ G (1) ( s, x, x ' ) , 0 ] , (6.11)  \nthe wave form nearly does not change with time, and it is exactly the Precursor. However, precisely the contribution from both the contour integral and G (2) ( s, x ' , x ) are nonzero, and this can be proven through numerical results (e.g.Fig.11), in which waveform slightly changes and even grows with time in Precursor.  \nWhen t > x ' + x , we arrive at the Ringdown stage, during which we have to choose the semicircle on the left-half plane for both terms, and now the result is  \nG ( t, x ' , x ) = Res [ G ( s, x, x ' ) , 0 ] + ∑ s n e s n t b ( s n ) ψ s n ( x ) ψ s n ( x ' ) 2 s n a ' ( s n ) + . . . , (6.12)  \nwhere s n = iω n and ω n is a certain QNM whose imaginary part is always positive, which means it contains contributions of all the QNMs. The ellipsis stands for those QNMs corresponding to higher-order poles or equivalently, Type-I instability. Now it's very clear that the evolution in this stage is exponential decay.  \nFor numerical calculation we follow the methods introduced in [35], that is, on the basis of the already attained L in (4.6), we can use the Hermite evolution matrix  \nA (4) = ( I -∆ τ 2 L + (∆ τ ) 2 12 L 2 ) -1 ( I + ∆ τ 2 L + (∆ τ ) 2 12 L 2 ) , (6.13)  \nto transform the field strength vector in the way of  \nu ( n +1) = A (4) u ( n ) , (6.14)  \nin which  \nu ( n ) = [ Φ( τ 0 + n ∆ τ ) Ψ( τ 0 + n ∆ τ ) ] . (6.15)  \nAlso, we set the similar initial data as in [22]:  \nΦ( σ, τ = 0) = exp( -( r ∗ ( σ ) -x 0 ) 2 8 r 2 + λ 2 ) , (6.16)  \n-∂ τ Φ | τ =0 = ∂ σ Φ | τ =0 , (6.17)  \nhere λ = 2 , x 0 ≫ 0 . Then we have the following results:  \n- 1. When p is larger the decay rate together with the oscillating rate is obviously much slower. Among the two important parameters, only p greatly contributes to the decay rate and the oscillation rate of the field strength, while q contributes just a little. This conclusion can be proven through comparing the two figures in Fig.11 and the comparison between the two figures in Fig.11 and the first figure in Fig.12. This conclusion warns that in real detecting the influence of quantum effect can easily be concealed by non-zero cosmological constant. Still larger q contributes to longer period in Precursor stage and slower decay rate, which is the most practically observable effect.\n- 2. The so called Forbidden Area mentioned in Fig.4 once again shows its important meaning in the last picture of Fig.12: in this area a non physical exponential explosion is unpreventable. Some modes with their imaginary parts negative is expected under this condition, which again illustrates the hyperboloidal method, at least through the construction introduced in this paper has lost effectiveness.\n- 3. The amplitude of the Ringdown stage almost has nothing to do with the specific position in x domain. This effect is particularly obvious when p is closer to 1. This shows a tendency towards average distribution of the total energy in the Precursor stage, and once this process is finished the whole evolution transfers to the Ringdown stage.  \nFigure 11: The time evolution of the Proca field on different positions under the parameter m = 0 , l = 1 , p = 0 . 99 and different q , with the initial data set in (6.16) and (6.17), in which x 0 = r ∗ ((10 r 0 + r + ) / 11) ≫ 0 .  \n<!-- image -->", 'VII. CONCLUSIONS AND DISCUSSION': "In this paper we study the quasinormal modes, their instability and pseudospectrum of the Proca field in quantum corrected qOS-dS spacetime. The qOS-dS solution is a solution under the LQG framework and have different characteristics to traditional RN-dS or Schwarzschild-dS solution. The correction can be summarized as two nonzero parameters, the relative ratio of the inner horizon radius to outer horizon radius q and the relative ratio of the outer horizon radius to cosmological horizon radius p . It's a meaningful question of how these two observably practical parameters relies on structural parameters α and Λ , and through  \nFigure 12: The time evolution of the Proca field on different positions under the parameter m = 0 , l = 1 , q = 0 . 5 and different p , with the initial data set in (6.16) and (6.17), in which x 0 = r ∗ ((100 r 0 + r + ) / 101) ≫ 0 .  \n<!-- image -->  \nparameterization this question now has an answer. After this was done we need to find the correct moving equation for Proca field in such a spacetime, and the axial perturbation's effective potential is obtained, which can be regarded as a small extension of Regge-Wheeler potential for vector fields. All these are our preparation work (Sec.II, Sec.III).  \nIn the first part of our formal work we construct the hyperboloidal framework in the classical way, that is, using the minimal gauge. We examine the effectiveness of this gauge under the most general parameter space and come to find that there really exists a certain parameter region where the spatial characteristic of the τ =constant hypersurface is violated. Following studies (Sec.VI) further illustrates the importance of guarantee the correctness of the construction. After the hyperboloidal framework and Chebyshev-Gauss-Labatto grid are successfully constructed, we present the QNMs under a certain parameter and certain numerical precision, and then confirm the astringency of the results by increasing the grid number together with increasing the machine precision.  \nIn the second part we study the instability of QNMs. We use two different methods to illustrate the instability. First is the parametric flow, in which we confirm that the QNMs all migrate in unrestricted lines when the correction from p or q becomes remarkable. When r + is fixed, no matter p or q , when they increase the average norms of QNMs will decrease. For parameter m , it hardly changes the imaginary part of QNMs while remarkably contributes to increase in the real part. This is quite different from previous studies on massive fields, in which m is believed to reduce the energy leak and lead to quasi-resonance. Actually m can even serve as a equivalent parameter as l . We then use the 2-norm of the matrix to define the pseudospectrum and reveal its characteristic through contour line plot, which show great accordance to the parametric migration flow.  \nIn the third part we study the evolution in the time domain both analytically and numerically. We first use the high-frequency approximation to obtain the expression of the two basic solution of the wave equation, and then the transfer coefficient and finally Green function. We illustrate how the contributions from different parts of contour integral influence the evolution. We then use numerical methods to clarify the impact caused by certain corrections. As expected, no power-law tail is found and p shows more importance than q . Still, quantum correction can influence the oscillating rate and the lasting period of the Precursor stage. These are all potential observables in related experiments in the future.  \nThere still remains plenty of questions to answer:  \n- 1. As m has been found to have nearly opposite influence in QNMs from that of quantum correction and cosmological constant, could it be possible that the massive term m precisely or approximately offsets the corrections caused by p or q ? If so, what is the specific parameter condition? This is very important in future detection if we want to determine the upper limit of parameters or oppositely, do parameter domain exclusion.\n- 2. The detailed properties of the gap between a very small p and p = 0 . In Fig.8, there are disturbing facts of immediate emerging new branches with rather small imaginary part when p approaches 0, which means the characteristic quasiresonance of massive fields in asymptotic flat spacetime are becoming remarkable gradually. Also, p → 0 and p = 0 majorly differs on whether there is branch cut and whether there is power-law tail, and this gap needs deeper understanding.\n- 3. How to attain the more precise approximation of the Green function as | s | → ∞ , as this determines how precisely we can predict the behavior in the Precursor stage.\n- 4. The strong cosmic censorship (SCC) hypothesis in qOS-dS spacetime has been tested under scalar perturbation in [18]. The conclusion there is actually based on the conclusion of [36] on scalar perturbation in RNdS, which means that is just an approximation. If we want to repeat this in Proca perturbation, we had better know more details about the relation between QNMs of vector field and SCC in quantum corrected spacetime.\n- 5. The entire work is based on semi-classical field theory, but strict conclusions still need to be done on complete LQG framework.", 'VIII. ACKNOWLEDGEMENTS': "We thank Li-Ming Cao for guidance, Liang-Bi Wu, Yu-Sen Zhou and Long-Yue Li for comments and very necessary assists in the technical instructions, Anıl Zengino˘glu for his fundamental works in hyperbolodital framework.  \n- [1] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116 , 061102 (2016), arXiv:1602.03837 [gr-qc] .\n- [2] E. H. T. Collaboration, The Astrophysical Journal Letters 875 , L6 (2019).\n- [3] P. Salucci and A. Burkert, The Astrophysical Journal 537 , L9 (2000).\n- [4] A. Burkert, The Astrophysical Journal 447 , L25 (1995).\n- [5] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Physics Reports 513 , 1 (2012), modified Gravity and Cosmology.\n- [6] E. Berti et al., Classical and Quantum Gravity 32 , 243001 (2015).\n- [7] E. Allys, P. Peter, and Y. Rodr'ıguez, Journal of Cosmology and Astroparticle Physics 2016 , 004 (2016).\n- [8] J. Beltr'an Jim'enez and L. Heisenberg, Physics Letters B 757 , 405 (2016).\n- [9] M. Minamitsuji, Phys. Rev. D 94 , 084039 (2016).\n- [10] H. Witek, V. Cardoso, A. Ishibashi, and U. Sperhake, Phys. Rev. D 87 , 043513 (2013).\n- [11] P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Phys. Rev. Lett. 109 , 131102 (2012).\n- [12] V. P. Frolov, P. Krtouˇs, D. Kubizˇn'ak, and J. E. Santos, Phys. Rev. Lett. 120 , 231103 (2018).\n- [13] M. Rahman and A. A. Sen, Phys. Rev. D 99 , 024052 (2019).\n- [14] M. Han, C. Rovelli, and F. Soltani, Phys. Rev. D 107 , 064011 (2023).\n- [15] T. Stachowiak and M. Szydłowski, Physics Letters B 646 , 209 (2007).\n- [16] J. Lewandowski, Y. Ma, J. Yang, and C. Zhang, Phys. Rev. Lett. 130 , 101501 (2023).\n- [17] C. Zhang, Y. Ma, and J. Yang, Phys. Rev. D 108 , 104004 (2023).\n- [18] C.-Y. Shao, C. Zhang, W. Zhang, and C.-G. Shao, Phys. Rev. D 109 , 064012 (2024).\n- [19] H. T. Cho, Phys. Rev. D 68 , 024003 (2003).\n- [20] A. Ishibashi and H. Kodama, Progress of Theoretical Physics 110 , 901 (2003).\n- [21] H.-P. Nollert, Phys. Rev. D 53 , 4397 (1996).\n- [22] R. G. Daghigh, M. D. Green, and J. C. Morey, Phys. Rev. D 101 , 104009 (2020).\n- [23] W.-L. Qian, K. Lin, C.-Y. Shao, B. Wang, and R.-H. Yue, Phys. Rev. D 103 , 024019 (2021).\n- [24] R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83 , 793 (2011).\n- [25] R. Konoplya and A. Zhidenko, Physics Letters B 609 , 377 (2005).\n- [26] Y. Hou, Z. Zhang, H. Yan, M. Guo, and B. Chen, Phys. Rev. D 106 , 064058 (2022).\n- [27] L.-M. Cao, J.-N. Chen, L.-B. Wu, L. Xie, and Y.-S. Zhou, (2024), arXiv:2401.09907 [gr-qc] .\n- [28] S. Chandrasekhar, 'Metric perturbation of ressiner-nodstorm solution,' in The Mathimatical Theory of Black Holes, edited by C.Wang (Higher Education Press, Beijing, 2018) pp. 216-225.\n- [29] J. L. Jaramillo, R. P. Macedo, and L. A. Sheikh, Phys. Rev. X 11 , 031003 (2021).\n- [30] Y. Yang, Z.-F. Mai, R.-Q. Yang, L. Shao, and E. Berti, 'Spectral instability of black holes: relating the frequency domain to the time domain,' (2024), arXiv:2407.20131 [gr-qc] .\n- [31] K. Destounis, R. P. Macedo, E. Berti, V. Cardoso, and J. L. Jaramillo, Phys. Rev. D 104 , 084091 (2021).\n- [32] S. Sarkar, M. Rahman, and S. Chakraborty, Phys. Rev. D 108 , 104002 (2023).\n- [33] R. F. Rosato, K. Destounis, and P. Pani, 'Ringdown stability: greybody factors as stable gravitational-wave observables,' (2024), arXiv:2406.01692 [gr-qc] .\n- [34] M. Casals and A. Ottewill, Phys. Rev. D 86 , 024021 (2012).\n- [35] M. F. O'Boyle, C. Markakis, L. J. G. Da Silva, R. Panosso Macedo, and J. A. V. Kroon, (2022), arXiv:2210.02550 [gr-qc] .\n- [36] P. Hintz and A. Vasy, Journal of Mathematical Physics 58 , 081509 (2017), https://pubs.aip.org/aip/jmp/articlepdf/doi/10.1063/1.4996575/13983995/081509 1 online.pdf ."}