Scalar Field Perturbation of Hairy Black Holes in EsGB theory

We investigate scalar field perturbations of the hairy black holes involved with spontaneous symmetry breaking of the global U(1) symmetry in Einstein-scalar-Gauss-Bonnet theory for asymptotically flat spacetimes. We consider the mechanism that black holes without hairs become unstable at the critical point of the coupling constant and undergo a phase transition to hairy black holes in the symmetry-broken phase driven by spontaneous symmetry breaking. This transition occurs near the black hole horizon due to the diminishing influence of the Gauss-Bonnet term at infinity. To examine such process, we introduce a scalar field perturbation on the newly formed background spacetime. We solve the linearized perturbation equation using Green's function method. We begin by solving the Green's function, incorporating the branch cut contribution. This allows us to analytically investigate the late-time behavior of the perturbation at both spatial and null infinity. We found that the late-time behavior only differs from the Schwarzschild black hole by a mass term. We then proceed to calculate the quasinormal modes (QNMs) numerically, which arise from the presence of poles in the Green's function. Our primary interest lies in utilizing QNMs to investigate the stability of the black hole solutions both the symmetric and symmetry-broken phases. Consistent with the prior study, our analysis shows that hairy black holes in the symmetric phase become unstable when the quadratic coupling constant exceeds a critical value for a fixed value of the quartic coupling constant. In contrast, hairy black holes in the symmetry-broken phase are always stable at the critical value. These numerical results provide strong evidence for a dynamical process that unstable black holes without hairs transition into stable hairy black holes in the symmetry-broken phase through the spontaneous symmetry breaking.


Introduction
Detecting gravitational waves from the merger of binary black holes first observed by the Laser Interferometer Gravitational-Wave Observatory (LIGO) [1] in 2015 was a significant physics breakthrough in recent decades.General relativity predicted the existence of gravitational waves, but their high frequency and weak nature made it challenging to confirm their existence experimentally, which took nearly a century.On the other hand, electromagnetic waves were observed just twenty years after their prediction from electromagnetic theory and have significantly impacted our lives.One of the primal objectives of LIGO is to test general relativity since it fails to account for dark matter, dark energy, and inflationary expansion [2][3][4].
Numerous alternative theories of gravity have been proposed to enhance general relativity [5][6][7][8][9].New terms or non-minimal coupling are introduced to modify the well-known Einstein-Hilbert Lagrangian.Models containing the Gauss-Bonnet term are of particular interest.
First, the Gauss-Bonnet term is a total derivative in four-dimensional spacetime [10].Second, although the gravity action contains higher curvature terms, it does not introduce higher derivative terms into the equation of motion, even when it is non-minimally coupled to a scalar field [5][6][7].Finally, the presence of the Gauss-Bonnet term is predicted from various string theory considerations [11,12].
This paper focuses on the Einstein-scalar-Gauss-Bonnet theory, where a scalar field couples to the Gauss-Bonnet term.In this theory, the evasion of the no-hair theorem was first studied in [13] (generalising results obtained in [14]) and soon after in [15] based on Bekenstein's argument [16,17].Paper [13] claimed that a positive coupling function is necessary for the evasion, but [15] showed that evasion is possible for both signs of the coupling function.These inconsistent results caused confusion and led to the conclusion that there is a privileged way to validate the evasion of the no-hair theorem.This would imply that the theorem loses its universal nature.The complete derivation for the evasion of the no-hair theorem was done in [18].It found that a surface term was omitted in the previous studies [13,15], resulting in inconsistent and incorrect results.The no-hair theorem is extended to spherical black holes in shift-symmetric scalar-tensor theory in [19].
Coinciding with the discovery of hairy black holes in ESGB theory [13], the mechanism of spontaneous scalarization was proposed to explain how black holes grow their scalar hair [20][21][22], which claims the tachyonic instability is responsible for this transition.However, it is found that the newly generated hairy black holes are dynamically unstable under the radial perturbation of scalar fields [23].This implies that this hairy black hole may not be the final stage in the evolution of a black hole.Alternative approaches were proposed later to guarantee the stability of the hairy black holes [24][25][26][27][28].
One promising approach in this regard utilizes spontaneous symmetry breaking (SSB) [29].The model contains a single complex scalar field with a global U(1) symmetry, which couples to the Gauss-Bonnet term: Here, G is the Gauss-Bonnet term, f is the coupling function, and V is commonly called the interaction potential.We only discuss the coupling function of the following form (1. 2) The interaction potential resembles the conventional scalar field potential in curved spacetime, with the Gauss-Bonnet term acting as a position-dependent magnitude.In this model, the Schwarzschild black hole becomes unstable under scalar field perturbation when the coupling constant reaches a critical value denoted as α Sch. .It was argued that this instability may lead to the formation of stable hairy black holes in symmetry-broken phases.Here the spontaneous symmetry breaking is involved in forming or evolving stable hairy black hole spacetimes from their counterparts without hairs.The Gauss-Bonnet term in the interaction potential peaks near the horizon and decreases rapidly as distance increases.Since the magnitude of the Gauss-Bonnet term is non-vanishing near the horizon, it is responsible for producing a nonvanishing expectation value for the scalar field that spontaneously breaks the U (1) symmetry.On the contrary, the contribution of the non-minimal coupling is negligible at infinity and so this process is irrelevant there.In the symmetry-broken phase, the equation of motion for Goldstone bosons is decoupled from other equations and due to the flux conservation only trivial solutions are physically acceptable.This application of spontaneous symmetry breaking to black hole systems is particularly intriguing, because it indicates that this mechanism, considered fundamental in our understanding of nature, might play a role in the formation and stability of hairy black holes.In the Standard Model, elementary particles gain mass through the Higgs mechanism and in condensed matter physics, spontaneous symmetry breaking is responsible for superconductivity.In cosmology, dark matter and dark energy may be composed of Goldstone bosons, created when symmetry is broken [30][31][32].
While the previous work in [29] identified a sufficient condition for (in)stability by analyzing the effective potential of the perturbed scalar field, this approach is not enough to prove the dynamical evolution of hairy black holes from their counterparts without hairs.To gain a more rigorous understanding of the stability for hairy black hole solutions, a detailed examination of their dynamical behavior is necessary.Calculating the quasinormal modes(QNMs) provides a powerful tool for such an investigation.
We aim to comprehensively analyze the scalar field perturbation on the fixed background of the hairy black hole solutions in the symmetric and symmetry-broken phases.We especially focus on calculating the QNMs to examine the instability of the hairy black holes.First, we solve the perturbed equation with Green's function.The function consists of three parts.G QNM describes the contribution of simple poles to scalar perturbations.G B describes the influence of the branch cut.G F contains all the other contributions coming from the part of the integration contour, which is two-quarter circles at infinity.We analytically obtained the solutions for G B in the low-frequency limit at spatial infinity (r → ∞) and null infinity (v → ∞) and demonstrated their solutions.As our background solutions are numerically generated, we utilized the numerical methods such as Chebysheve polynomials to compute the QNMs to calculate G QNM .We proved that the hairy black holes in the symmetric phase become unstable at a certain value of α (denoted as α crit.), yielding the positive value for the imaginary part of the frequency, while the hairy black holes in the symmetry-broken phase are stable for all valid values of α, yielding the negative value for the imaginary part of the frequency.As discussed in [29], we also found α Sch.≈ α crit. .These numerical results suggest that black holes without hairs undergo a transition to hairy black holes in the symmetrybroken phase.
This paper is organized as follows.We introduce the gravity model and hairy black hole solutions in Section 2. Section 3 delves into the Green's function for scalar field perturbations on the hairy black hole background.We discuss the full solution structure of the scalar field perturbation, demonstrate some properties of QNMs, and show the late time behavior by calculating the branch cut in Green's function.Next, we discuss the numerical method employed to calculate QNMs in Section 4. Sections 5 and 6 focus on calculating the QNMs for hairy black holes.Section 5 demonstrates the hairy black holes in the symmetric phase and examines their QNMs (Section 5.3).Section 6 follows the same arguments for the hairy black holes in the symmetry-broken phase.Finally, Section 7 presents a summary of the key findings.

Hairy Black Holes by Spontaneous Symmetry Breaking
We consider the Einstein-scalar-Gauss-Bonnet theory in four-dimensional asymptotically flat spacetime as follows: where κ 2 = 8πG the following scalar field coupling function is employed where φ(r) is a complex field and α and λ are coupling constants.We assume that α can have real values, but λ only takes positive ones.This Lagrangian respects the global U (1) symmetry where χ is a constant.We are specially interested in forming stable hairy black holes from black holes without hairs, suggested in [29] where the transition mechanism is realized by the spontaneous symmetry breaking.This will happen if the potential of the scalar field is invariant under a certain symmetry and possesses local minima.We treat the scalar field coupling with the Gauss-Bonnet term as an interacting potential V = −f (φ * , φ)G.Depending on the values of the parameters in f , the interacting potential can have either a single-well or double-well shape at every point in the spacetime.This allows us to investigate hairy black holes in two distinct phases: symmetric and symmetry-broken phases.These phases are determined by the sign of α or the vacuum value of the interacting potential.The symmetric phase occurs when the scalar field near the horizon resides at either the "global" minimum (when α < 0) or the "local" maximum (when α > 0) of the interacting potential.In contrast, the symmetry-broken phase is characterized by the scalar field near the horizon settling at the "global" minimum (α > 0) of the interacting potential.In later sections, we will construct hairy black hole solutions for both the symmetric and symmetry-broken phases.
We use the spherically symmetric and diagonal form of the metric ansatz The equations of motion take the following form Here and it is expanded as (2.7) A complex scalar field can be decomposed into two real scalar fields as follows : Thus, we replace the complex field with the real scalars.We only consider the case of asymptotically flat spacetime.Consequently, near infinity, the metric behaves as where i, j = 1, 2. The coefficient A 1 = −2M is identified with the black hole mass M .The scalar charge is defined as follows 12) The definition immediately results in φ i,1 = D i .For simplicity, we consider φ 1 = φ 2 hereafter.

Green's function for Scalar Field Perturbation
The hairy black hole solutions satisfying the equations of motion (2.5) -(2.6) are hard to obtain analytically.However, they can be established with numerical methods.On these hairy black hole background, presented in Section 5 and 6, we consider scalar field perturbations and understand their solutions using Green's function.

Scalar Field Perturbations
Let us consider a perturbation of the scalar field described by (2.6) The equation takes the following form where the subscript of f indicates the variation with respect to those variables.We decompose the complex scalar with (2.8) and employ the tortoise coordinate dx = 1 The perturbation equation becomes To solve the equation for δφ 1 (t, x), we impose boundary conditions reflecting the classical nature of black holes: an incoming condition at the horizon (x → −∞) and an outgoing condition at infinity (x → ∞).These conditions imply that a black hole in asymptotically flat spacetime acts as an open system, and energy carried by δφ 1 (t, x) dissipates.Consequently, the system is non-Hermitian, meaning its energy eigenvalues are complex rather than real.In Hermitian systems, energy is conserved and energy eigenvalues are real, leading to normal modes forming a complete set.In contrast, an open system has complex energy eigenvalues, and the corresponding solutions yield QNMs.They do not form a complete set and experience exponential decay or divergence over time depending on the imaginary part of the frequency.
To search for solutions of (3.4) one can use the Laplace transformation on δ φ1 (s, x) for real s and positive t.
If δφ 1 (t, x) is bounded, one can analytically continue δ φ1 (s, x) into the complex half-plane Re(s) > 0. After the transformation, the perturbation equation takes the following form The source term J (s, x) is determined by the initial conditions at t = 0.The equation is an inhomogeneous ordinary second-order differential equation that can be solved using the Green's function.By definition, the Green's function is a solution of the following equation: where x ′ and x are the distances to the source and the observer from the horizon, respectively.
Once the Green's function is obtained and the initial conditions are given, one obtains the following solution of the perturbation equation Here the Green's function is constructed by the solutions of the homogenous equations φ(±) with J (s, x) = 0 from (3.7) as follows Then the Green's function takes the form where W x (s) indicates the Wronskian computed in the x-coordinate.
To investigate the general properties of δ φ1 (s, x), we temporarily neglect the contribution of the initial condition to the source term J (s, x ′ ) and instead, focus on the solution of the Green's function.We obtain the Green's function in t-plane by performing an inverse Laplace transformation To calculate this integral, one shall construct the closed contour integral in the complex splane.Since t ≥ 0, the integral diverges for Re(s) > 0, which is the right half of the s-plane.Thus, we make a closed contour integral for negative Re(s), which is the left half of the splane.One should also consider a branch cut for real negative s, so the integration contour should not cross the branch cut.Consequently, the Green's function consists of three parts: where ) The first term G QNM (t, x, x ′ ) describes the contribution of poles within the contour γ R , the second term G B (t, x, x ′ ) arise from the branch cut, and the last term G F (t, x, x ′ ) is originated from the two quarter-circles at infinity.When a scalar perturbation is introduced to the background spacetime, it generates an initial wave that propagates both inward toward the black hole and outward toward a distant observer.If the initial perturbation carries significantly higher energy compared to the black hole's effective potential, denoted by V eff (r) in (3.5), the outward propagating wave can effectively ignore the influence of this potential.Consequently, this high-energy component traveling outward reaches the observer at a time approximately equal to t ≃ x − x ′ .This solution is described by the term G F (t, x, x ′ ).However, while this high-frequency signal might offer insights into the early response of the black hole to the perturbation, calculating the exact behavior of the wave interacting with the potential is highly complex and remains an open problem.
At t ≃ x + x ′ , the part of the low-frequency scalar perturbation traveling towards the black hole bounces off the effective potential V eff (r) near the horizon and propagates back to the observer.This reflected signal is described by the term G QNM (t, x, x ′ ), and their basic properties are discussed in Section 3.2.There are several methods known to compute the QNMs [33][34][35][36][37][38].Because only numerical background solutions are known, we use numerical methods to compute their QNMs.We discuss this in detail in Section 4. We present the numerical results of these QNMs in Sections 5 and 6, enabling a rigorous investigation of stability for our black hole solutions.
At late time t ≫ x + x ′ , the QNMs gradually fade away, and only the decaying modes of the signal persist.This late-time behavior arises from the outgoing wave that traveled to a large distance and then scattered back towards the observer due to the spacetime curvature at large distances.Consequently, the behavior at these late times depends on the asymptotic form of the effective potential in (3.5) at infinity.This signal can be seen in G B (t, x, x ′ ), and we analytically calculate this quantity for our hairy black hole solutions in section 3.3.

Some Properties of Quasinormal Modes
In this subsection, we closely investigate the physical properties of the QNMs in (3.16) which is the contour integral that is only contributed by the poles.We assume that the observer is far away from the source (x > x ′ ), and so our Green's function in (3.12) is written as Since φ+ and φ− are analytic within the contour γ R , poles can only appear because of the Wronskian zeros.In other words, the function has poles in points where W x = 0. We suppose that W x only has a set of simple poles (s = s n ) and expand W x near these poles The residue theorem yields both Green's function and the scalar field perturbations to take the sums over all poles ) where and At the same time, the condition W x = 0 in (3.13) also lead its integration to yield where c 0 is a function depending only on s.This implies that φ+ and φ− are dependent.Another distinct characteristic of QNMs is their non-normalizability under the Klein-Gordon inner product: Expanding (3.7) near infinity and changing s = −iω, we obtain where x ∼ r from (3.3) is used, the source term is neglected since the variable change to ω indicates the Fourier transformation.Separating variables as δ φ1 (ω, r) = e −iωt δ φ1 (r), the asymptotic solution reads where c 1 and c 2 are integration constants, the second line is expanded for large r.In this solution, terms e −iω(t−r) r −iωA 1 = e −iω(t−r+A 1 log r) and e −iω(t+r) r iωA 1 = e −iω(t+r−A 1 log r) describe an outgoing and incoming waves respectively.We require the outgoing boundary condition at infinity, so we set c 2 = 0. Since the frequency is discretized, the solution is written as Plugging this to the Klein-Gordon inner product, we obtain where Σ t is a hypersurface with a timelike normal vector n t .Since we are working with asymptotic solutions, we employed dΣ t = dΣ n t = −g tt g rr r 2 sin θdrdθdϕ ≃ r 2 sin θdrdθdϕ by substituting the asymptotic expansions from (2.9)- (2.11).When n = m, the inner product becomes where we made a large r expansion at the second line after integration.As mentioned above, A 1 is the black hole mass, with 2M = −A 1 > 0. For stable black holes, the imaginary part of the quasinormal mode frequency, denoted by Im(ω n ), is negative (i.e., Im(ω n ) < 0).Consequently, the Klein-Gordon inner product diverges at the boundary of infinity when evaluated for QNMs.This divergence implies that QNMs are non-normalizable and do not form a complete set.

Branch Cut
Here, we compute the part of the Green's function (3.17) that accounts for the branch cut, following the approach in [39].This term takes the following form: Let us return to the r-coordinate and introduce a new variable ψ(s, r) In this variable the equation (3.11) reads: Expanding (3.34) to the large r regime by using (2.9)-(2.11), it yields Here, l 1 accounts for both the angular momentum, coming from the spherical harmonic contribution Y m l and the leading order contribution from the metric and the scalar coupling function: Whittaker functions M and W solve the equation(3.35) and they have properties of In turn, the Wronskians take the following form Using these properties, the Green's function in the s-space is computed as where the Wronskian is Consequently, the branch cut contribution to the Green's function reads (3.45) In terms of the r coordinate, this expression becomes where we used the Wronskian computed in the r coordinate Since we consider the low-frequency regime (s ∼ iω for Fourier modes), we take the limit A 1 s → 0 and φ 1,1 s → 0. This leads the Green's function to be where 1 F 1 (2rs) = 1 F 1 (l + 1; 2l + 2; 2rs) is the confluent hypergeometric function of the first kind.Firstly, we consider the low-frequency limit rs ≪ 1 and r ′ s ≪ 1.In the limit t → ∞ with r fixed, we obtain Our late-time solution aligns with the leading order term of the Schwarzschild solution at large r in [39,40], when we substitute −A 1 for 2M .
In the limit v → ∞ (t → ∞ with t/r fixed), we can take r → ∞ and r ′ s ≪ 1.Consequently, we obtain the following expression where v = t + x ∼ t+r and u = t − x ∼ t − r and v → ∞ is taken at the last step.This result also agrees with the radiative case for no initial static field at null infinity in [40] by replacing the mass of hairy black holes.
-12 -Since the background black hole solutions are numerical, we use a numerical method, namely, the Chebyshev spectral method, to calculate the QNMs.This subsection discusses the employed numerical methods initially introduced in [38].We present results for black holes in the symmetric and the symmetry-broken phases in Sections 5 and 6 correspondingly.One needs to set up a grid to realize a space to use numerical calculations.Here we employ Chebyshev polynomials T n (y) defined by which introduces unequispacing grids [41,42].If one sets cos θ = y, then the polynomials become The Chebyshev polynomials form a complete basis, so any function defined on an interval [−1, 1] admits the Chebyshev expansion: The polynomials are also orthogonal with respect to the following inner product ⟨f, g⟩ implying Due to the limit of the numerical capability, we approximate the Chebyshev polynomial T n (y) to N -th degree polynomial We generate (N + 1) Chebyshev points y n = cos( nπ N ), n = 0, 1, • • • , N accordingly.The derivative of the function of f (y) is expressed as where M is a matrix that is made up with the coefficients from the derivative of the Chebyshev polynomial with respect to y on the basis of the Chebyshev polynomial T n (y).It takes a form of  The higher derivatives of f (y) are obtained by For example, if the perturbed equation has a form of we can express the functions as follows We rearrange the above expression as Here I is the (N + 1) × (N + 1) identity matrix and The subscript of β indicates the power of ω.Then we redefine α, β 1 , and β 2 to take a form of where For the case of equations containing higher orders of ω, we can generalize α and β as -14 - where n max is the highest power of ω.By using the command "Eigenvalues[{α, β}]" in Mathematica, we can compute the quasinormal frequencies ω.We here employ the Chebyshev zero points which generate N Chebyshev points [42].Then we obtain N values of ω.Because of the numerical errors we repeat the numerical calculation with different numbers of grids and changing the accuracy and precision.If these calculations produce the same values of ω within the tolerance of 10 −2 , we consider that those values of ω are reliable candidates.(In our numerical calculation, the tolerance ranges from 10 −35 to 10 −2 .)Then, we check the eigenfunctions for the selected ω in the previous step.If they satisfy our boundary conditions, we consider those values of ω to be correct.
We first describe our numerical solutions on the grids based on these properties and numerical techniques.Since the Chebyshev polynomials are defined for −1 ≲ y ≲ 1, we transform the coordinate r, which ranges from the black hole horizon to infinity, into a finite interval (z h ≤ z ≤ z ∞ ) using the substitution r = 1/z.Subsequently, we discretize z into z n by rescaling where z ∞ = 0.Then, the Chebyshev points are shifted to 0 ≲ z n ≲ 1.
To incorporate the incoming and outgoing boundary conditions into our numerical calculations, we use ingoing Eddington-Finkelstein (EF) coordinates.These coordinates have a distinct advantage: when solving the perturbation equation near the boundaries (both horizon and infinity), the asymptotic solutions naturally separate into incoming and outgoing modes.As we will see in more detail later, the incoming modes are generally normalizable and exhibit smooth behavior.Conversely, the outgoing modes are non-normalizable and oscillate rapidly as they approach the horizon or infinity.Because of such rapid oscillations, the wavelength of the outgoing mode becomes much smaller than the grid spacing.Hence, our numerical calculations will not be able to capture it effectively.This limitation is beneficial near the horizon because we aim to impose the incoming boundary condition there.However, it becomes problematic at infinity, where we require the outgoing boundary condition that involves the non-normalizable and rapidly oscillating modes.To make the outgoing modes numerically treatable, we can redefine the outgoing mode to become normalizable.This transformation will automatically render the incoming mode non-normalizable and rapidly oscillating.By strategically redefining the field in this manner, we can effectively impose the boundary conditions for both incoming and outgoing QNMs on our numerical calculations.

Quasinormal Modes for Hairy Black Holes in the Symmetric Phase
We numerically generate hairy black hole solutions in the symmetric phase for both positive and negative values of α.Based on prior theoretical expectations [29], hairy black holes beyond a certain positive value of α are anticipated to be unstable.To verify this, we calculate their QNMs and demonstrate that they indeed become unstable under scalar field perturbations at a critical value of α, denoted as α crit. .

Background solution
In order to use our numerical method, we transform it to the ingoing EF coordinate as follows where z = 1/r and v is defined as Then, the Gauss-Bonnet term is written as where ′ is the derivative with respect to z. Equations of motions (2.5)-(2.6)are expressed in terms of φ 1 and φ 2 as follows where (5.8) Asymptotic flatness is required at infinity.Under this condition, the metric and scalar fields are expanded as follows: where all coefficients are constant.As before, we identify the coefficient A 1 as ADM mass of the black hole such that A 1 = −2M , and φ i,1 is the scalar charge, φ i,1 = Q i .Hereafter we set φ 1 = φ 2 for simplicity.
For a regular black hole to exist, specific boundary conditions must be met near the horizon: where ϵ = z h − z is the expansion parameter and A h , B h , φ ih and φ ih,1 (i = 1, 2) are constants.We also set κ 2 = 1/2.The following relations between these constants ensure the regularity of the metric and scalar fields: As seen above, the near horizon expansion of the fields is determined uniquely by two boundary values, which are A h and φ 1h = φ 2h = φ h .We discovered that φ ′′ (r h ) diverges when the value inside the root becomes zero.To avoid this, we impose the regularity conditions: (5.15) Thus, this regularity condition constrains the values of parameters (α, λ) and the scalar field value near the horizon (φ h ).In Figure 1, we present the phase space for these parameters at a fixed value of λ.The shaded regions enclosed by the solid lines represent valid parameters for hairy black hole solutions.Points lying on the solid lines themselves are not considered valid.To obtain the solutions, we impose the boundary value of φ h for given values of α and λ and numerically solve the equations of motion (5.4)-(5.7)from the near horizon (z = z h − ϵ) to infinity (z = ϵ).Here, A h is determined by the boundary conditions at infinity.To perform numerical calculation, we fixed z h = 1 and ϵ = 10 −50 and used "NDSolve" implemented in Mathematica with precision up to 300.High precision and accuracy in constructing the background solution are crucial for reliable computation of quasinormal frequencies.This is because numerical errors in the background solution directly propagate to the calculated frequencies.As mentioned in [29], numerical solutions were not found in the entire valid parameter space but found beyond the certain value of α, which is denoted as α sol. in Figure 1.The absence of solutions in the region α < α sol. is discussed in [29].We present the numerically generated hairy black hole solutions in the symmetric phase for two different scalar field values at the horizon (φ h ).Given the value of λ = 0.1, Figure 2 shows the numerical solutions for the metric and scalar field with φ h = 0.01, comparing with the Schwarzschil black holes.To highlight the contrasting behavior of the solutions, Figure 3 zooms in on a smaller region of the z-axis.This magnified view clearly reveals a significantly more pronounced anisotropy (A/B ̸ = 1) in the metric for large negative α values compared to positive α values.

Scalar Field Perturbation
We consider a linearized perturbation of the complex scalar field around the background solution so the linearized equation becomes (5.17This deviation, though subtle, leads to noticeably distinct behavior in their QNMs.Conversely, for α = 3.5, the hairy black hole solutions show a much weaker deviations from the Schwarzschild metric across the radial distance.
Employing the following ansatz for scalar field perturbation and using the linearized scalar field equation (5.17) is written as and (5.21) Let us investigate the asymptotic solution of (5.20) at both boundaries.Near the horizon, the perturbation equation is expanded as and plugging Φ(z) ∼ (z h − z) p = ϵ p into above it expands For the equation not to be singular, the first term should vanish.This requires p to take values: which yield δφ(v, z) = e −iωv const. (5.25) . (5.26) -20 - The first solution corresponds to the incoming mode.The second one is written as δφ(v, z) = e −iω(v−2z * +2c) = e −iω(u+2c) , (5.27) where the outgoing null vector u = t − z * is used and z * is expanded near the horizon such as where c indicates the first three terms in the second equality.This indicates that the second solution in (5.24) corresponds to the outgoing mode.Also, as shown in (5.26), the phase of this outgoing solution diverges and oscillates very fast when z approaches z h in (v, z)-coordinate.
Due to the limitations of numerical calculations on a grid, the rapid oscillations of the outgoing modes are not captured.Consequently, we can only impose the incoming boundary condition at the horizon, effectively excluding the outgoing modes from the numerical calculation.Near infinity the perturbation equation (5.20) is expanded and the solution takes the following form To rigorously check this solution, we generalize the form to be .31)and expand (5.29) near z ∼ 0. Then we obtain (5.32) Taking z → 0 limit, the first and second line require p = 1 + 2iA 1 ω and p = 1 respectively and they yield The second term having c 2 describes the incoming mode and is normalisable.The first term having c 1 is rewritten as where u = t − z * is an outgoing null vector and z * in (5.2) is expanded in z → 0 limit This indicates that the first term describes the outgoing mode but is non-normalizable and oscillates more rapidly as z approaches 0 in (v, z)-coordinate.Here, we want to impose the outgoing boundary condition.This situation is problematic.To resolve this we redefine Φ(z) as and the new field Φ(z) becomes Then, the incoming mode becomes non-normalizable and is not captured by the numerical calculation due to the fast oscillation, while the outgoing mode becomes normalizable and behaves smoothly.Rewriting the perturbation equation (5.20) in terms of Φ(z), it yields By adopting the numerical method explained in Section 3.4, we will numerically calculate the quasinormal frequency ω in (5.38) and display their values in the following subsection.

Our Numerical Results
We first consider the Schwarzschild black hole, which is the solution with φ = 0, and examine its stability according to the value of α.By computing the quasinormal frequency, we found that a positive imaginary part of the frequency, Im[ω], emerges at a critical value, α Sch.≈ 0.1815 as shown in Figure 4.This indicates that the corresponding system becomes unstable α ≥ α Sch. .Here, we want to argue that this instability triggers a transition from the Schwarzschild black hole to a hairy black hole, and the hairy black hole will be formed not in the symmetric phase but in the symmetry-broken phase.To support this argument, we show that the hairy black holes with α ≥ α Sch. are unstable in the symmetric phase but stable in the symmetry-broken phase.This section focuses solely on the quasinormal frequencies of hairy black holes in the symmetric phase for various α values.The stability analysis of hairy black holes in the symmetry-broken phase will be addressed in Section 6.
To compute the eigenvalues, we performed calculations with grids of size 100 and 140, resulting in 100 and 140 eigenvalues, respectively.To identify reliable eigenvalues, we considered only those for which both the real and imaginary parts agreed at least to the second decimal places (i.e.tolerance of 10 −2 ) between the two grid sizes.Finally, we plotted the corresponding eigenfunctions associated with these converged eigenvalues to verify -22 - that they satisfy the boundary conditions.Figure 5 shows the quasinormal frequencies for α = −0.2,−0.8, −0.9, −1, −1.1, −1.2, −1.3 and −1.4,generated on the background solution with φ h = 0.01 for a fixed value of λ = 0.1.Figure 6 demonstrates the quasinormal frequencies for α = 0.12, 0.1815, 0.
.4 and 2.8 on the same background.The detailed numerical values for the QNMs presented in Figure 5 and 6 can be found in Appendix A.
Our results show that the imaginary parts of the QNM frequencies are all negative for the negative values of α.In Figure 5, the QNMs exhibit a well arranged line of spectrum at α = −0.2.As α decreases towards more negative values, a progressive shift is observed in the higher-order QNMs.The fifth QNM exhibits this behavior first, moving inwards and upwards relative to the main QNMs' array.Subsequently, at larger negative α values, the sixth QNM starts displaying the same behavior.This sequence of inward and upward shifts with increasing negative α forms the additional array of the QNM spectrum inside the outer array of the QNMs.However, pushing α to even more negative values leads to a depletion of higher-frequency QNMs, resulting in a sparser spectrum at α = −1.4. Figure 3 reveals that the metric exhibits increasing anisotropy (A/B ̸ = 1) with more negative values of α.This deviation from isotropy (A/B = 1) is suspected to be a contributing factor to the emergence of the exotic (cases of α = −0.8 to −1.2) and diminishing (cases of α = −1.3 and −1.4) QNM spectrum observed at large negative α values.
For the positive values of α, we found that the imaginary parts of the QNM frequencies become positive at α = 0.1815 (denoted as α crit.).This positive imaginary part indicates that the hairy black hole becomes dynamically unstable in the regime α ≥ α crit. .We demonstrated the corresponding eigenfunctions to ensure these eigenvalues and examined their behaviors in Figure 14 in Appendix B. To investigate the angular dependence of the QNM spectrum -23 -in this case, we further checked the QNM behavior for different values of l.The results are presented in Figure 15 in Appendix C. Figure 6 illustrates that as the value of α increases, lower-frequency QNMs (with smaller imaginary parts) transition to having positive imaginary parts.As the value of α increases, our current numerical approach struggles to capture high-frequency QNMs with negative imaginary value.These high-frequency QNMs might become sparser, meaning they are fewer in number and more challenging to detect with our current methods.Also the nature of these modes could be affected by the metric's instability, causing them to lose their characteristic of oscillating and decaying over time and potentially transitioning into unstable modes themselves.For α = 3.5, the QNMs exhibit purely positive imaginary frequencies.This implies that the oscillations of the perturbed scalar fields vanish and are only damped or growing at this value.We also computed the QNMs frequencies for α = −0.3,0.184, 0.1848 and 0.4 for φ h = 0.1 and λ = 0.1 and plotted them in Figure 7. Table 2 presents the numerically computed QNM frequencies.Black colored numbers indicate a high degree of accuracy, with a numerical error less than 10 −5 between the results obtained from grids of size 100 and 140.Conversely, gray colored numbers signify potential numerical inaccuracy in those digits.In this case, we found that the critical value of α appears at α crit.= 0.1848.Since these critical values of α are involved in dynamical stability, their values differ from those computed in [29], which are sufficient conditions.

Quasinormal Modes for Hairy Black Holes in the Symmetry-Broken Phase
This section demonstrates the numerically generated hairy black holes' solution in the symmetrybroken phase and shows their QNMs.

Background solution
When α is positive, the interaction potential V forms degenerate vacuums, and the stable minima are described by where the vacuum states are labeled by the parameter β.These ground states do not respect the symmetry of the Lagrangian, which indicates that symmetry is spontaneously broken in the vacuum.We expand a field around a ground state v by reparameterizing it as follows Here, the new variables σ(z) and θ(z) are physical fields because they are excitations above the vacuum.Then we rewrite the Lagrangian in terms of σ(z) and θ(z)  -25 -  -26 -  For the scalar field flux to be finite, c 2 should be zero [29].
As we consider the regular black hole solutions, we impose the boundary condition near the horizon as follows where the near horizon expansion of the metric and the scalar field σ(z) take the following forms from the equations of motion The two independent variables A h and σ h and the coupling constants α and λ fully describe the boundary conditions.The regularity condition requires which limits the values of α, λ and σ h .Figure 8 illustrates the valid parameter space for hairy black hole solutions at fixed values of λ.The shaded regions bounded by the solid lines represent valid values of α and σ h , excluding points that fall exactly on the solid lines.Within the parameter values in the phase space, we numerically generate hairy black hole solutions in the symmetry-broken phase with the same precision and accuracy used in the symmetric phase.For the choice of α = 0.01 with the fixed value of λ = 0.1, we plot the solutions in Figure 9 and Figure 10.By increasing the value of α, we observed that the numerical solutions become increasingly difficult to find for values of α exceeding α sol , as discussed in [29].

Scalar Field Perturbation
We consider a linearized perturbation of the real scalar field as follows σ → σ(z) + δσ(v, z, θ, ϕ).(6.17) Applying c 2 = 0 in (6.12) to the scalar field equations (6.6), the linearized part becomes Using the following ansatz and (5.19), the linearized equation is written as where The equation structure of (6.21) remains unchanged from (5.20).Thus we utilize (5.38), substituting fφ * φ with f σσ , to test the instability of the QNMs.

Our Numerical Results
We employed the same strategy used in the symmetric phase to ensure the convergence of the eigenvalues.Calculations were performed on grids of 100 and 140, resulting in 100 and 140 eigenvalues, respectively.To identify reliable eigenvalues, we considered only those for which both the real and imaginary parts agreed within a tolerance of 10 −2 between the two grid sizes.Here, we only used hairy black hole solutions with a fixed value of λ = 0.1.To illustrate the behavior of quasinormal frequencies for varying α values, we first computed them for a background solution with σ h = 0.01.Representative results are presented for α = 0.01, 0.2, 0.34, and 0.38 in Figure 11, and for α = 0.42, 0.46, 0.50, 0.54, 0.58, 0.60, 0.62, and 0.66 in Figure 12.The numerical values for Figure 12 are presented in Appendix A. The QNM spectrum in Figure 12 for increasing positive values of α exhibits behaviors similar to those observed in Figure 5 for increasing negative values of α.Different from the symmetric phase with φ h = 0.01, the imaginary parts of quasinormal frequencies are all negative, which indicates the stability of hairy black holes.We also examined the quasinormal frequencies for σ h = 0.1 and plotted them for various values of α (α = 0.01, 0.1, 0.2, and 0.22) in Figure 13, which shows all negative values for the imaginary parts of quasinormal frequencies.We further investigated the quasinormal frequencies with the background solutions with σ h = −0.01 and σ h = −0.1 for various values of α.It turns out that they are stable for the linear scalar perturbation and the quasinormal frequencies.
-32 -  -33 - -34 - We have argued that black holes without hairs grow their scalar hairs by spontaneous symmetry breaking under global U(1) symmetry in Einstein-scalar-Gauss-Bonnet theory, following the suggestion in [29].We investigated the scalar field perturbation on the hairy black hole spacetimes to provide evidence for this process.We have fully analyzed the solution for the linearized scalar field equation by using Green's function, which consists of three parts (G F , G B , and G QNM ).Our main focus was the calculation of QNMs (G QNM ), which provides a powerful tool to study the stability of concerning systems.The branch cut contribution (G B ) differs from the Schwarzschild black hole case only by the mass term.Since our hairy black hole solutions are numerically generated, we employed the numerical methods to compute the QNMs.
We firstly considered the Schwarzschild black hole, which is the solution without scalar fields, and examined the QNMs of scalar perturbation on this background.We found a critical point at which the imaginary part of the quasinormal frequency becomes positive, indicating the system's instability.Then, we numerically constructed hairy black hole solutions in the symmetric and symmetry-broken phases.We assumed the scalar field always lies near the vacuum in each phase, so their near horizon values are small.Our calculations of the QNMs show that the hairy black holes in the symmetric phase also have a positive imaginary value of the quasinormal frequency at nearly the same critical point as the Schwarzschild black hole.However, the hairy black holes in the symmetry-broken phase always produce a negative imaginary value of the quasinormal frequency in the parameter ranges that we are concerned with.This numerical result indicates that the Schwarzschild black holes undergo a phase transition to the ones in the symmetry-broken phase by spontaneous symmetry breaking.Thus, their final stage through evolution becomes stable.
Our investigation focused on the scalar field perturbation within the probe limit, neglecting the potential backreaction on the metric.Since the background scalar field has a non-trivial contribution to the spacetime geometry, the perturbation itself will induce a corresponding change in that geometry.This backreaction effect can potentially influence the stability of the system.Thus, it is crucial to consider this effect in future work.Furthermore, studying the full non-linear time evolution of this process would be much more intriguing, but this work presents a significantly greater technical challenge.6.
spectrum for l = 0 to 5 and their numerical data are presented in Figure 15 and Table 8 respectively.

Figure 1 :
Figure 1: Phase space of α and φ h for λ = 0.1 in the symmetric phase )

( a )Figure 2 : 5 Figure 3 :
Figure 2: Hairy black hole solutions with φ h = 0.01 for various values of α.The black dotted line denotes the solutions of Schwarzschild black hole.The small regimes with the callout index (a)∼(f) are plotted in Figure 3.

Figure 8 :
Figure 8: Phase space for α and σ h with fixed values of λ in the symmetry-broken phase

( a )Figure 9 : 6 Figure 10 :
Figure 9: Hairy black hole solutions with σ h = 0.01 for various values of α.The black dotted line denotes the solutions of Schwarzschild black hole.The small regimes with the callout index (a)∼(f) are plotted in Figure 10.

Figure 14 :
Figure 14: The eigenfunctions corresponding to the quasinormal frequencies of the hairy black hole with φ h = 0.01 at α = 0.1815 and λ = 0.1.The subscript of Φ(z) shares the same number of the QNM data presented in Table6.

Table 1 :
QNMs of Schwarzschild black holes for serval values α