Perturbations of massless external fields in Horndeski hairy black hole

In this paper, we study the propagations of external fields in Horndeski theory, including the scalar field, electromagnetic field and Dirac field. We extensively explore the quasinormal frequencies, time evolution, greybody factors and emission rates of those massless perturbing fields by solving the corresponding master equations in the Horndeski hairy black hole. With the use of both numerical and analytical methods, we disclose the competitive/promotional influences of the Horndeski hair, spin and quantum momentum number of the external fields on those phenomenal physics. Our results show that the Horndeski hairy black hole is stable under those perturbations. Moreover, a larger Horndeski hair could enhance the intensity of energy emission rate for Hawking radiation of various particles, indicating that comparing to the Schwarzschild black hole, the Horndeski hariy black hole could have longer or shorter lifetime depending on the sign of the Horndeski hair.

Recent observational progress on gravitational waves [1][2][3] and black hole shadow [4,5] further demonstrates the great success of Einstein's general relativity (GR).Yet it is unabated that GR should be generalized, and in the generalized theories extra fields or higher curvature terms are always involved in the action [6][7][8].Numerous modified gravitational theories were proposed, which indeed provide a richer framework and significantly help us further understand GR as well as our Universe.Among them, the scalar-tensor theories, which contain a scalar field ϕ as well as a metric tensor g µν , are known as the simplest nontrivial extensions of GR [9].One of the most famous four-dimensional scalar-tensor theory is the Horndeski gravity proposed in 1974 [10], which contains higher derivatives of ϕ and g µν and is free of Ostrogradski instabilities because it possesses at most second-order differential field equations.Various observational constraints or bounds on Horndeski theories have been explored in [11][12][13][14][15][16].
Horndeski gravity attracts lots of attention in the cosmological and astrophysical communities because it has significant consequences in describing the accelerated expansion and other interesting features, please see [17] for review.Moreover, Horndeski theory is important to test the no-hair theorem, because it has diffeomorphism invariance and second-order field equations, which are similar to GR.In fact, hairy black holes in Horndeski gravity have been widely constructed and analyzed, including the radially dependent scalar field [18][19][20][21][22][23][24][25][26][27][28] and the time-dependent scalar field [29][30][31][32][33][34].However, the hairy solution with scalar hair in linear time dependence was found to be unstable, and so this type of hairy solution was ruled out in Horndeski gravity [35].Later in [36] , the no-hair theorem was demonstrated not be hold when a Galileon field is coupled to gravity, but the static spherical black hole only admits trivial Galileon profiles.Then, inspired by [36], the authors of [37] further examined the no-hair theorem in Horndeski theories and beyond.They demonstrated that shift-symmetric Horndeski theory and beyond allow for static and asymptotically flat black holes with a nontrivial static scalar field, and the action they considered is dubbed quartic Horndeski gravity where χ = −∂ µ ϕ∂ µ ϕ/2 is the canonical kinetic term, Q i (i = 2, 3, 4, 5) are arbitrary functions of χ and Q i,χ ≡ ∂Q i /∂χ, R is the Ricci scalar and G µν is the Einstein tensor.In particular, very recently a static hairy black hole in a specific quartic Horndeski theory, saying that Q 5 in the above action vanishes, has been constructed in [38] Here, M and Q are the parameters related to the black hole mass and Horndeski hair.The metric reduces to Schwarzschild case as Q → 0, and it is asymptotically flat.From the metric (2), it is not difficult to induce that for arbitrary Q, r = 0 is an intrinsic singularity as the curvature scalar is singular, and f (r) = 0 always admits a solution r + = 2M which indicates a horizon at r + = 2M .In addition, when Q > 0, r + = 2M is the unique root of f (r) = 0 , so it has a unique horizon, i.e., the event horizon for the hairy black hole as for the Schwarzschild black hole.While for −2M < Q < 0, f (r) = 0 has two roots: one is r + = 2M indicating the event horizon, and the other r − denotes the Cauchy horizon which is smaller than the event horizon.r − increases as Q decreases, and finally approaches r + as Q → −2M , meaning the extremal case.This hairy black hole and its rotating counterpart attract plenty of attentions.Some theoretical and observational investigations have been carried out, for examples, the strong gravitational lensing [39], thermodynamic and weak gravitational lensing [40,41], shadow constraint from EHT observation [42], superradiant energy extraction [43] and photon rings in the black hole image [44].However, the propagation of external field in the Horndeski hairy black hole is still missing.To fill this gap, here we shall explore the quasinormal frequencies (QNFs), time evolution, greybody factors and emission rates by analyzing the massless external fields perturbations (including the scalar field, electromagnetic field and Dirac field) around the Horndeski hairy black hole.
QNFs of a field perturbation on the black hole are infinite discrete spectrum of complex frequencies, of which the real part determines the oscillation timescale of the quasinormal modes (QNMs), while the complex part determines their exponential decaying timescale.They dominate the signal of the gravitational waves at the ringdown stage and are one of the most important characteristics of black hole geometry.The interest of QNMs in more fundamental physics can be referred to the reviews [45][46][47].Mathematically, QNFs depend solely on the basic three parameters of the black holes, i.e., the mass, charge, and angular momentum.However, if there are any additional parameters that describe the black hole, such as the hairy parameter in this framework, those parameters will also have prints on the QNMs spectrum.Even though the propagations of external field in a black hole background seems less related to the gravitational wave signals, they still might provide us with important insights about the properties of the Horndeski hairy black holes, such as their stability and the possible probe of the characterized parameters of black holes.This is our motivation to investigate the external fields QNMs of the hairy black hole solution (2).The goal is to study the influences of the hairy parameter Q on the QNFs signature for the massless scalar field, electromagnetic field and Dirac field perturbations, respectively.To this end, we will use both the WKB method and the matrix method to numerically obtain the QNFs, and also exhibit the time evolution of the perturbations in the time domain.
The other goal of this work is to study the impact of Horndeski hair on the energy emission rate of the particles with spin = 0, 1 and 1/2, respectively, and the greybody factors of their Hawking radiation from the Horndeski hairy black hole.The grey-body factor measures the modification of the pure black body spectrum and it is equal to the transmission probability of an outgoing wave radiated from the black hole event horizon to the asymptotic region [48].It significantly describes information about the near-horizon structure of black holes [49].So, one can evaluate the energy emission rate of Hawking radiation with the use of the greybody factor [50].It is known that the Hawking radiation spectrum and its greybody factor are very sensitive to the modifications of GR.So they at least provide important sources of physical consequences for the modifications in the formulation of the black hole.Thus, we could expect that the Horndeski hair will leave prints on the Hawking radiation spectrums as well as the greybody factors.It is noted that the study of the energy emission rates and greybody factors of the particles with spin= 0, 1 and 1/2 requires one to solve the master equations of the scalar, electromagnetic and Dirac perturbing fields on the hairy black hole background.This process is similar to that we do when calculating the quasinormal modes, only with different boundary conditions: the latter requires a purely outgoing wave at infinity and a purely ingoing wave at the event horizon, while the former allows ingoing waves at infinity.This paper is organized as follows.In section II, we show the master equations for the test massless scalar, electromagnetic and Dirac fields in the Horndeski hairy black hole, and analyze the properties of their effective potentials.In section III, we calculate the QNM frequencies of the perturbing fields with both the WKB method and matrix method, and then match the behaviors of the perturbations in the time domain.In section IV, by solving the corresponding master equations, we evaluate the greybody factor and the energy emission of Hawking radiation for various particles.The last section contributes to our conclusions and discussion.Throughout the paper, we will set c = ̵ h = G = 1.Moreover, in a convenient way, we will fix M = 1/2 and denote Q/M → Q in all the computations.

II. MASTER EQUATIONS AND EFFECTIVE POTENTIALS OF THE PERTURBING EXTERNAL FIELDS
In this section, we will show the master equations of various massless external fields, including the scalar field, electromagnetic field and Dirac field around the Horndeski hairy black hole.The influences of Horndeski hair Q and angular quantum number ℓ on the corresponding effective potentials of the perturbations will be analysed.

A. Scalar field perturbation
The propagation of massless scalar field Φ in the Horndeski hairy black hole satisfies the Klein-Gordon equation where g is the determinant of the black hole metric (2).By taking the ansatz where ω is the frequency of scalar field perturbation and m is the azimuthal number, we can separate (3) and obtain the radial master equation Here, the prime denotes a derivative w.r.t.r, and the effective potential is where ℓ = 0, 1, 2, ⋯ is the angular quantum number.Under the tortoise coordinate the equation ( 5) can be written into Schrodinger-like form The behavior of effective potential V sc (r) with various samples of Q and ℓ are shown in FIG. 1.It can be obviously seen that, for each case with different Q and ℓ, the potential functions are always positive outside the event horizon.
The absence of a negative potential well may give a hint that the black hole could remain stable under the massless scalar field perturbation.These plots also show that the effective potential always has a barrier near the horizon, which is enhanced by increasing the values of Q and ℓ.FIG.1: The effective potential Vsc(r) for the massless scalar field perturbation.In the left plot we fix ℓ = 0 and tune Q, while we fix Q = −0.5 and tune ℓ in the right plot, respectively.

B. Electromagnetic perturbation
The propagation of electromagnetic field in the Horndeski hairy black hole background satisfies the Maxwell equation where is the field strength tensor, and A µ is the vector potential.In order to separate the Maxwell equation, we take A µ as the Regge-Wheeler-Zerilli decomposition [51,52] where Y ℓm = Y ℓm (θ, φ) are the scalar spherical harmonics with the angular quantum number, ℓ, and azimuthal number, m, respectively.By substituting Eq.( 10) into Eq.(9), we will obtain two decoupled radial equations which can be uniformed into the master equation with for axial modes with odd parity (−1) ℓ+1 , [iωh ℓm (r) + dj ℓm (r) dr ] for polar modes with even parity (−1) ℓ . ( Again under the tortoise coordinate ( 7), the Schrodinger-like form of Eq.( 11) reads as where the effective potential is The potential function for the electromagnetic field perturbation is depicted in FIG. 2 which shows that both larger Q and ℓ give higher potential barrier, similar to the case in scalar field perturbation.
FIG. 2: The effective potential of electromagnetic perturbation V EM (r).We fix ℓ = 1 and tune Q in the left plot, while we fix Q = −0.5 and tune ℓ in the right plot, respectively.

C. Dirac field perturbation
The propagation of massless fermionic field in the Horndeski hairy black hole background is governed by In this equation, (ω µν ) a is the 1-form spin connections; (e µ ) a is a rigid tetrad defined by (e µ ) a = √ g µν (dx µ ) a and its dual form is (e µ ) a = (e ν ) b g ab η µν with η µν the Minkowski metric; Γ a is the curved spacetime gamma matrices, which connects the flat spacetime gamma via Γ a = (e µ ) a Γ µ .After working out the spin connections of the metric (2), we expand the Dirac equations as To proceed, we choose the representation of the flat spacetime gamma matrices as [53] where σ i are the Pauli matrices.Then considering the Dirac field decomposition [53] where the spinor angular harmonics are we can obtain two radial master equations where κ ± = ∓(j + 1 2 ) for j = l ± 1/2.The Schrodinger-like equations under the tortoise coordinate take the forms where the effective potentials are It was addressed in [54] that the behaviors of V I Dirac and V II Dirac usually are qualitatively similar because they are super-symmetric partners derived from the same super potential, so one can choose one of them to proceed without any loss of generality.Thus, in the following study, we will concentrate on the master equation ( 22) with V I Dirac , which is plotted in FIG. 3 for some references of Q and ℓ.
Comparing these three effective potentials in FIG.1-FIG.3 for various perturbing external fields, we can extract the following properties.(i) The Horndeski hair has the same effect on the maximum value of various potentials, i.e., a larger Q corresponds to a larger and narrower barrier, from which we expect the same influence on the QNFs of those perturbations.(ii) The effect of ℓ on the behavior of various potentials in Horndeski hairy black hole are similar to that in the Schwarzschild black hole [55].(iii) The barrier of the effective potential for perturbing field with higher spin seems larger and wider.We could expect that the features of effective potentials could be reflected in the QNFs.FIG.3: The effective potential of Dirac perturbation V I Dirac (r).We fix ℓ = 0 and tune Q in the left plot, while we fix Q = −0.5 and tune ℓ in the right plot.

III. QUASI-NORMAL MODE FREQUENCIES OF VARIOUS PERTURBATIONS
In this section, we shall compute the quasi-normal frequencies of the Horndeski hairy black hole under the massless scalar, EM and Dirac fields perturbations by solving the master equations ( 8), ( 13) and ( 22) with the boundary conditions: ingoing wave (∼ e −iωr * ) at the horizon and the outgoing wave (∼ e iωr * ) at infinity.We will employ both the WKB method and matrix method to guarantee that we focus on the fundamental mode with the node n = 0, and also to pledge the credibility of our results.Moreover, in order to directly analyze the essence of QNFs in the propagation of various perturbations, we will also study the time evolution of the perturbation fields with the use of time domain integration.Instead of the tedious details in the main part, we briefly review the main steps of WKB method, matrix method and domain integration method of our framework in appendixes A-C, because all the methods are widely used in the related studies.Especially, we use the Padé series P ñ m [56] in WKB method to improve and reconstruct the WKB correction terms and transform it into a continued-fraction-like form.Here ñ and m are the order numbers of the numerator and denominator series (see [57]).
A. Q− dependence Firstly, we analyze the influence of the Horndeski hair parameter Q on the QNM frequencies of the lowest ℓ modes in various external field perturbations.The results are listed in TABLE I (also depicted in FIG. 4) in which we also calculate the relative error between two methods defined by The QNFs obtained from the matrix and WKB-Padé methods agree well with each other.For various perturbations, the imaginary part of QNFs, Im(ω), keeps increasing as Q decreases, and it is always negative even in the extremal case with Q = −1.It means that the Horndeski hairy black hole is dynamically stable under those external fields perturbation with the lowest-lying ℓ.Moreover, by comparing the QNFs for various perturbations, we find that for the field with larger spin, the Im(ω) is larger (see also the left plot of FIG. 4 ).Thus, the perturbation field with a higher spin could live longer than the one with a lower spin because the damping time τ d for a wave field is related with the QNF by τ d ∼ 1/| − Im(ω)|.However, the real part of QNFs, Re(ω), for all the perturbations decreases as Q decreases, meaning that smaller Q suppresses the oscillation of the perturbations.Similar to the imaginary part, Re(ω) is larger for a perturbing field with higher spin.The effects of Horndeski hair and the spin of fields on the QNFs can be explained by their influences on the corresponding effective potentials as we described in the previous section.

B. ℓ− dependence
Now we move on to study the effect of the angular quantum number for various perturbations.To this end, we focus on Q = −0.5.The results are shown in TABLE II and FIG. 6.Following are our observations: (i) For small ℓ, it has a relatively strong influence on the Re(ω), while the influence on the Im(ω) is weak.We observe that as ℓ increases, Im(ω) for the EM field tends to slightly decrease, while the Im(ω) for scalar and Dirac fields slightly increase.This means that the effect of growing ℓ will shorten the lifetime of the EM perturbation, but it can instead extend the lifetime of the scalar and Dirac perturbations.Similar phenomena can also be observed in the Schwarzschild black hole [58], but here we find that the introducing of the Horndeski Q has print on the effect of ℓ.In detail, positive Q will enhance the effect while negative Q suppresses this effect.These findings in QNFs can be verified from FIG. 7 and its comparison to FIG. 5. (ii) As ℓ increases, the results from WKB methods and matrix match better and better, and the relative error becomes smaller than 10 −6 .This is because the WKB-Padé method is essentially a semi-analytic approximation which works better for an analytical form with higher ℓ [59].(iii) When we further increase ℓ, the gap among the QNFs for all the massless perturbing fields tends to be smaller.This is because for large ℓ, the dominant terms in all the effective potentials ( 6), ( 14) and ( 24) are the terms ∝ ℓ 2 , which have the same formula in all cases.And finally the gap will vanish in the eikonal limit ℓ ≫ 1 which we will study in the next subsection.The fundamental (n = 0) QNMs of various massless perturbation modes with different angular momentums obtained by WKB method and matrix method and their relative errors.Here we fix Q = −0.5.

C. QNFs in eikonal limit
Cardoso et al proposed that in the eikonal limit (ℓ ≫ 1), the real part of QNM frequency for a static spherical black hole is connected with the angular velocity of the circular null geodesics while the imaginary part is connected with the Lyapunov exponent [60].Then the real part of the QNMs in the eikonal limit was further related to the shadow radius of a static black hole as ω Re = lim ℓ≫1 ℓ R sh [61], and more recently this connection was extended into the rotating black holes [62].This correspondence may originate from the fact that the perturbing waves could be treated as massless particles propagating along the last timelike unstable orbit out to infinity, but deeper research deserves to be done for further understanding.
In this subsection, we compare the quasinormal spectrum obtained from the master equations with the spectrum computed directly from the geometric-optics approximation formula, which is given by [60] with , and for the Horndeski hairy black hole (2).Ω c is the angular velocity of a massless particle geodesically moving on a circular null orbit with radius r = r c , and λ LE is the Lyapunov exponent, where the radius r c is given by the positive root of the equation 2f c = r c f ′ c .In TABLE III, we list the QNFs obtained from wave analysis and the geometric-optics approximation for fixed ℓ = 40.QNFs for various perturbing fields converge to be the value obtained from geometricoptics correspondence, because in the limit ℓ ≫ 1, all the perturbed wave equations reduce to the analytical equation describing the geodesic motion of the massless particle in the Horndeski hairy spacetime.Then in order to check the effect of Q, we plot Ω c and λ LE as functions of Q in FIG. 8. We see that both of them grow monotonically as Q increases, indicating a smaller imaginary part but a larger real part of QNFs for ℓ ≫ 1.The effect of Q on the QNFs is already reflected for small ℓ ∼ 1 as we disclosed in the previous subsections.On the other hand, the results in previous subsection show that the increasing Q can enhance the contribution of the spin so that QNFs of perturbing fields with different spins emerge a larger bifurcation (see FIG. 4 ).However, in the eikonal limit the perturbing fields with different spins tend to possess a same QNFs, which indicates that the contribution of field spin is diluting as ℓ increases.In order to analytically understand the balance effect of Q and ℓ on the spin contribution on the QNFs, we apply the Newman-Penrose formalism to construct general spherical symmetric Teukolsky equations for an arbitrary field spin s in static spherical-symmetric metric with ds 2 = −f (r)dt 2 + f −1 (r)dr 2 + r 2 (dθ 2 + sin 2 θdφ 2 ), which separate into angular and radial parts as with ∆ ≡ r 2 f (r), m the azimuthal number, ϵ s and A sℓ listed in TABLE IV.For more details on the formula derivation, readers can refer to [63][64][65].It is straightforward to check that the above equations can reduce to the Teukolsky equation in static limit derived in [63].To proceed, we reform R s = ∆ −s/2 Ψ s /r and work in tortoise coordinate dr * = (r 2 /∆)dr , thus the radial equation ( 30) is transformed into the wave-like equation V sℓ (r) = i s ω r 2 d dr ( It is noted that the Teukolsky radial equations are different from the corresponding master equations ( 8), ( 13) and (22) when setting s = 0, 1, 1/2, because they are not in the form of canonical wave equations.However, it was addressed in [66] that the Teukolsky equations can be brought into the corresponding master equations under certain transformation.The QNFs are significantly determined by the effective potential and dominated by the last term in V sℓ .In the Horndeski hairy black hole (2), the influence from the spin can be amplified through the coupling s 2 Q 2 originating from the term s 2 ∆ ′2 /4r 4 , which eventually causes the bifurcation or, saying the wider "fine structure" in the quasinormal spectrum.While in the eikonal limit ℓ >> 1, the term ∆A sℓ /r 4 → ∆ℓ 2 /r 2 shall become dominant in the potential, so that the effect of spin will be suppressed.This subsequently causes the degeneracy for both the null circular orbits and the quasinormal spectrum of massless fields with different spins.

IV. GREYBODY FACTOR AND HAWKING RADIATION
It is known from Hawking's paper [50] that a particle with negative energy measured by an observer at infinity can physically exist inside the black hole since there the Killing vector is spacelike.Thus, during the pair production at the vicinity of the horizon, the particle with positive energy can escape to the observer and leave the negative one to fall into the singularity.This effect causes the so-called Hawking radiation, whose power spectrum is shown to be a literally black-body spectrum.However, due to the existence of the potential barrier outside the black hole, the Hawking radiation is not totally transparent for the observer at infinity, so what the observer detects in fact is a grey-body spectrum because the particles could be scattered by the potential barrier.
In order to quantize the scattering process for various particles, one should first calculate the transmission coefficient, which is also defined as the greybody factor.Here we will solve the master equations ( 8), ( 13) and ( 22) by the scattering boundary conditions which permit the ingoing wave at infinity differing from the case of calculating QNFs.The boundary conditions that can equivalently describe the scattering process for a particle emitted from the horizon read as where T ℓ and R ℓ are denoted as the transmission and reflection coefficients for the angular momentum ℓ mode, satisfying Since from section II we know that each effective potential exhibits a potential barrier decreasing monotonically towards both boundaries, thereby, we can again employ the WKB method described in appendix A to calculate the coefficients.Subsequently, we have the greyfactor |A ℓ | 2 defined as [59,67] where K can be obtained from the WKB formula Here, V (r 0 ) denotes the maximal potential locating at r 0 for various spins, the second derivative in V is with respect to the tortoise coordinate, and Λ i are the higher order WKB correction terms.It is noted that the WKB formula for determining the grey-body factors is well known to provide reasonable accuracy for further estimating the energy rate of Hawking radiation.However, this approach may not be suitable for the cases with very small ω, which imply almost complete wave reflection with negligible contributions to the total energy emission rate, so we set a cutoff for ω in the numeric to make our results reliable.
With the above preparation of methodology, we depict the grey-body factor as a function of ω with different Horndeski hair Q for various fields in FIG. 9. Two remarkable properties we can extract from the figure.(i) Similar as in Schwarzschild black hole (Q = 0), the increasing of spin s and orbital angular momentum ℓ could make the curves of the grey-body factor shift towards a larger ω in Horndeski hairy black hole.This means that for the emission particle with larger ℓ and s, the lowest frequency with which the particle penetrates the potential barrier would raise up such that the barrier would tend to shield more low-frequency particles for the observer at infinity.(ii) Regarding to the effect of the hairy charge Q on |A ℓ | 2 , we observe that the increasing of |Q| would not make the curve of |A ℓ | 2 shift distinctly, instead, it makes the curve branch out between two fixed frequencies.It means that the Horndeski hair would not change the lowest frequency for a particle transmitting the potential barrier.Moreover, by comparing the dotted, solid and dashed curves with the same color for all the fields, it is obvious that for larger Q, the greybody factors is smaller for fixed ω, which means that more number of particles is reflected by the corresponding effective potential.This observation is reasonable because for larger Q, all the effective potentials barrier is higher (see FIG. 1-FIG.3), so the particles are more difficult to penetrate.With the greybody factor in hands, we can further estimate the energy emission rate of the Hawking radiation of the Horndeski hariy black hole via [50] where ± in the denominator denote the fermions and the bosons created in the vicinity of the event horizon, the Hawking temperature T H is and the multiples N ℓ is known as Note that the formula (36) only works when the system can be described by the canonical ensemble, which is fulfilled under the assumption that the temperature of black hole does not change between two particles emitted in succession [50].
The energy emission rate (EER) for scalar, Dirac and Maxwell fields as a function of frequency with samples of ℓ and Q are shown in Fig. 10.We can extract the following features.(i) The general behavior of EER is that as ω grows, the EER first increases till it reaches a maximum at ω = ω p where ω p depends on the parameters, and then decays to be zero.This is because in the right side of (36), there exists a competition between the greybody factor in the numerator and the exponential term in the denominator.Though FIG. 9 shows that |A ℓ | 2 grows as ω increases, the growth of |A ℓ | 2 is only more influential for ω < ω p , while for ω > ω p the growing of e ω/T H plays the dominant role and makes the EER decrease.And further increasing ω, |A ℓ | 2 tends to the unit, but the denominator increases exponentially, which causes the EER to decay exponentially.(ii) In each plot, we see that the EER for various fields in the black hole with a larger Q is stronger.The result is reasonable because for fixed ω, as Q increases, the greybody factor becomes smaller and the exponential term in the denominator also decreases due to the growth of Hawking temperature (37).The joint contributions result in the stronger intensity of EER for larger Q. (iii) By comparing the two plots for each field, we observe that EER for particles with higher orbital angular momentum would be suppressed in both Schwarzschild and Honrdeski hairy black holes.In addition, this suppression effect is more significant for the particles with a higher spin.
Finally, we numerically integrate the EER over ω and obtain the total EER, dE/dt, for various fields as the function of Q in FIG.11.For the static hairy black hole in the current Horndeski gravity, the positive Horndeski hair Q will enhance the total EER around the black hole, which may lead to a higher speed of evaporation and a shorter lifetime for this kind of black hole according to the discussion in [68].Therefore, in terms of observation, this intuitively would mean that, such hairy small or even medium black hole with a large positive Q in the early universe may have disappeared due to the high evaporation rate.While for negative Q, one would expect that the Horndeski hairy black hole has a lower evaporation rate.Especially, in the extremal case with Q = −1, the total EER approaches zero as expected, meaning that such extremal black holes almost do not evaporate and thus they are often considered to live forever if they are in complete isolation, similar to the case in extremal RN black hole [69].

V. CONCLUSION AND DISCUSSION
In this paper, we investigated the quasinormal frequencies and Hawking radiations of the Horndeski hairy black hole by analyzing the perturbations of massless fields with spins 0 (scalar field), 1/2 (Dirac field) and 1 (electromagnetic field), respectively.The starting points of both aspects are the master equations of the perturbing fields.For the quasinormal mode spectral analysis, we employed three methods: WKB method, matrix method and time domain integration; while in the Hawking radiation part, we used the 6−th order WKB method to determine the greybody factor.
Our results show that under the massless perturbations of the scalar field, Dirac field and electromagnetic field, the Horndeski hairy black holes are dynamically stable in terms of the frequency and time domain of QNMs.Similar as in Schwarzschild black hole [55], the massless field with higher spin has a larger imaginary part of quasinormal frequency, so the related perturbation can live longer in the hairy black hole background.This effect can be enhanced by having a stronger negative Horndeski hair.Moreover, the real part of QNFs for all the perturbations increases as Q increases, meaning that the larger Horndeski hair enhances the oscillation of the perturbations.In addition, in Horndeski hairy black hole, the mode with a larger ℓ has a shorter lifetime for electromagnetic field perturbation, but it can survive longer for the scalar and Dirac perturbations, which are similar to those in Schwarzschild black hole.But this effect of ℓ on the QNFs would be enhanced by a positive Q while suppressed by a negative Q. Nevertheless, for large enough ℓ till the eikonal limit (ℓ ≫ 1), the QNFs for various perturbations tend to be almost the same value, of which the imaginary (real) part decreases (increases) as Q increases.Focusing on the general Teukolsky equations for arbitrary spin in the current hairy black hole, we explained the balance effects of the Horndeski hair, the spin and the quantum angular momentum on the QNFs in an analytical way.
Our studies on Hawking radiation shew that the intensity of energy emission rate for various fields is stronger for a larger Horndeski hair, which indicates that the Horndeski hair could remarkably influence the evaporation rate of black hole.Thus, we argue that such black hole with large positive Horndeski hairy charge (if exists) in the early universe may have disappeared due to the strong radiation, while comparing to black hole in GR, Horndeski hairy black hole with negative hair could have a lower evaporation rate.
It would be more practical to extend our study into the gravitational perturbations, i.e., the massless spin-2 field, and we believe that it deserves an individual consideration due to the difficulty in reducing the master wavelike equation in the current theory.However, our findings from the test external fields could shed some insights on the related physics on gravitational perturbations.For example, the QNFs of gravitational fields perturbation in eikonal limit could be the same as our findings because the behavior of the QNM spectrum for external and gravitational fields is usually known to be qualitatively the same, independent of the spin of the field in this limit.In addition, our findings about the influence of the field's spin on the QNFs and the energy emission rate could also provide a good reference in this scenario.
Moreover, considering that the overtone modes may play important roles in the GW [70,71], it would be interesting to further study the QNFs for the overtone modes.Additionally, the extension to the external fields with mass could be another interesting direction.At least the propagation of the massive scalar field was found to behave differently from the massless scalar field [72,73].We hope to perform these studies in the near future.
To proceed, we consider the coordinate transformation to bring the integration domain from r ∈ [r + , ∞] to x ∈ (0, 1].Further by implementing the function transformation we can reform the equation (41) into A 2 (x)χ ′′ (x) + A 1 (x, ω)χ ′ (x) + A 0 (x, ω)χ(x) = 0, (47) where the expressions of A i (i = 0, 1, 2) are straightforward and will not present here.Subsequently, the boundary conditions are simplified as With the reformed master equation ( 47) and the boundary conditions (48) in hands, we can follow the standard steps of matrix method directly to obtain the eigenvalue ω in it.Following is the main principles of the matrix method [75][76][77].One first interpolates N grid points x 1 = 0 < x 2 < x 3 < ⋅ ⋅ ⋅ < x N = 1 in the interval x ∈ [0, 1], then by carrying out Taylor expansion of χ at anyone of these grid points x k (k = 1, 2, ⋯, N ), one has Finally, by taking x = x j (j = 1, 2, ⋯, k − 1, k + 1, ⋯, N ), one can get a matrix equation where ∆F = (χ(x 1 ) − χ(x k ), χ(x 2 ) − χ(x k ), ⋯, χ(x k−1 ) − χ(x k ), χ(x k+1 ) − χ(x k ), ⋯, χ(x N ) − χ(x k )) It is more convenient to express χ ′ (x k ) and χ ′′ (x k ) as the linear combination of the function value of χ at each grid point using the Cramer rule, where the matrix M i (i = 1, 2) is constructed with replacing the i'th column of matrix M by ∆F .In this way, one can finally transform the master equation ( 47) into a matrix equation M (ω)F = 0, (52) equations and effective potentials of the perturbing external fields A. Scalar field perturbation B. Electromagnetic perturbation C. Dirac field perturbation III.Quasi-normal mode frequencies of various perturbations A. Q− dependence B. ℓ− dependence C. QNFs in eikonal limit IV.Greybody factor and Hawking radiation V. Conclusion and discussion I. INTRODUCTION

FIG. 4 :
FIG.4: Quasi-normal frequencies as a function of the hairy charge at the low-lying angular quantum number, i.e., ℓ = 0 for scalar & Dirac fields and ℓ = 1 for EM field.

FIG. 5 :
FIG.5: Time evolution of the lowest-lying mode for the perturbing scalar field with spin s = 0 (left), Dirac field with s = 1/2 (middle) and EM field with s = 1 (right), respectively.

1 , ℓ = 2 FIG. 10 :
FIG.10: Energy emission rate as a function of frequency for various fields with spin s = 0 (top), s = 1/2 (middle), s = 1 (bottom) around Horndeski hairy black hole.The left column describes the fields with the corresponding lowest ℓ while the right column is for their second lowest ℓ.

TABLE III :
QNFs with ℓ = 40 for various perturbing fields obtained via different methods.
FIG.8:The angular velocity and the Lyapunov exponent as functions of hairy charge Q with fixing ℓ = 40.

TABLE IV :
Variables A Sℓ & ϵs dependent on the field spin in spherical symmetric Teukolsky equations.