{'fraction_non_alphanumeric': 0.04695434123915996, 'fraction_numerical': 0.013039389823424929, 'mean_word_length': 4.210801393728223, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Schwarzschild scalar wigs: spectral analysis and late time behavior

Using the Green's function representation technique, the late time behavior of localized scalar field distributions on Schwarzschild spacetimes is studied. Assuming arbitrary initial data we perform a spectral analysis, computing the amplitude of each excited quasi-bound mode without the necessity of performing dynamical evolutions. The resulting superposition of modes is compared with a traditional numerical evolution with excellent agreement; therefore, we have an efficient way to determine final black hole wigs. The astrophysical relevance of the quasi-bound modes is discussed in the context of scalar field dark matter models and the axiverse.

On the other hand, black holes are a natural consequence of general relativity. They appeared for the first time as a purely mathematical result, but they were later recognized as real astrophysical objects. Presently, we believe black holes represent the final fate of sufficiently massive stars, and supermassive black holes probably reside at the center of most galaxies [29].
Since both scalar fields and black holes seem to inhabit our Universe, a natural question arises: How do scalar fields and black holes interact with each other? Static, spherically symmetric, asymptotically flat black holes with a non-trivial scalar field distribution cannot exist in nature [30,31]; however, these results do not rule out extremely long lived transient phenomena [32,33].
The massive Klein-Gordon equation on a black hole background has been studied extensively.
For a Schwarzschild black hole Ternov et al. [34] found the exis-tence of quasi-bound modes (also called quasi-stationary modes or quasi-resonances), describing scalar field configurations surrounding the black hole for some period of time; see also Refs. [35,36] for later work on these states. Similar modes were found on a Kerr black hole background [34,37,38]. In the rotating case, there is a specific range of parameters for which these modes are exponentially growing in time, implying that a massive scalar field on a rotating Kerr black hole is unstable. See Refs. [39][40][41][42] for rigorous results on this instability, and [43][44][45] for recent numerical studies. Based on Leaver's continuous fraction method [46], Konoplya and Zhidenko [47] computed the quasi-normal frequencies of a massive scalar field on a Schwarzschild spacetime background, and mentioned conditions under which quasi-bound states could arise. For a generalization to the Kerr-Newman case, see Ref. [48]. Finally, we mention that Dolan and collaborators [49,50] generalized the computation of quasi-bound modes to different types of fields, vectorial or fermionic.
Motivated by the possibility of describing the dark matter in the galactic halo by the coherent excitation of an ultra-light scalar field [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28], including the supermasive black hole [51][52][53][54], in previous work we analyzed the time scale for the quasi-bound states on a supermassive, nonrotating black hole background [32,33]; see [55] for a recent generalization to the Dirac field. In particular, we showed in Ref. [33] that these states can actually survive for cosmological times, as long as the product of the black hole and the scalar field mass is sufficiently small. Furthermore, we provided strong evidence in Ref. [33] that these quasi-bound modes are generic. This conclusion was reached by numerically evolving arbitrary initial data and Fourier-transforming the solution with respect to time. The spectrum revealed a clear excitation of the quasi-bound frequencies.
In this paper we use the Green's function representation technique [56] to determine the late time behavior of localized scalar field distributions on a Schwarzschild background. Within this technique, the Green's function of the problem is decomposed into three contributions, corresponding to different contours in the complex frequency plane. The contribution we are focusing on in this work is the one arising from the residua of the poles of the Green's function, which describe the excitations of the quasi-bound states. This allows us to predict, given arbitrary initial data for the scalar field on a spacelike hypersurface, the amplitude of each quasi-bound state in the solution by computing a simple, one-dimensional integral. Furthermore, we show by comparison with results obtained from numerical evolutions that for large times (within the time scales of our simulations) the behavior of the solution is well approximated by the corresponding superposition of quasi-bound sates; therefore, the remaining two contributions to the Green's function become negligible for those times. Note that this situation is different than the case of a massless scalar field propagating on a Schwarzschild background, where the oscillating part of the solution, described by a superposition of quasi-normal modes, is soon taken over by the slower, polynomial tail decay. This difference resides in the fact that the quasi-bound modes have a decay rate which is much smaller than the quasi-normal modes, so that they dominate the solution for much larger time.
At this point we would like to stress that our definition of "large" or "late" times in this article is the following: we consider time scales which are very large compared to the Schwarzschild radius of the black hole, but which are within the lifetime of the quasi-bound states. In particular, we do not claim that for t → ∞ the solution is described by a simple superposition of quasi-bound states since in this limit the contribution from the tail part is expected to play an important role. However, for the physical scenario we are interested in, it is enough to consider large times which, according to our definition, could be larger than the age of the universe as we explained above.
Our results provide the basis for a generalization to the case of a rotating black hole. Several lines of work have been developed in the context of scalar and other fields on Kerr background; there are those which consider unstable modes in the vicinity of the Kerr black hole [57,58], and even possible radiation emission due to such instabilities [59,60]. The continuation of our work, however, will focus on a somewhat different direction, namely the scalar field in the Kerr background in connection with dark matter halos.
The remaining part of this work is organized as follows. In Sec. II we briefly review the Green's function method and the corresponding decomposition of the solution in three different contributions. Next, in Sec. III we discuss the case of a potential formed by a delta and a step function, as a toy model for which the modes can be represented analytically. In particular, we determine the quasi-bound frequencies in the limit for which the height of the step is very small compared to the amplitude of the delta barrier and show that they are located very close to the resonant band, resulting in very small decay rates. Then, we compute the excitations of each quasi-bound state, compare the resulting superposition to the results obtained from a numerical evolution, and find a very good agreement at late times.
In Sec. IV we apply the Green's function technique to the case of a massive scalar field on a Schwarzschild background. Here the modes and quasi-bound frequencies are computed semianalytically via a matching algorithm from which the excitation amplitudes can be determined. In Sec. V we compare the semianalytic results with those obtained from numerical evolutions of the Cauchy problem, and as in the case of the toy model, for late times we find an excellent agreement between the two approaches. Then, in Sec. VI we use our semianalytic solution to extract some physical properties of the quasi-bound states that might be relevant in the astrophysical context, and finally in Sec. VII we present our conclusions.
We also include two appendixes. In the first one we generalize the Green's function technique to time foliations different than the one arising from standard Schwarzschild coordinates. In the second appendix we provide some estimates on the location of the quasibound frequencies for the toy model example and the time after which the solution can be described by a superposition of the corresponding modes.

II. GREEN'S FUNCTION REPRESENTATION AND SPECTROSCOPY
In this section we briefly review the Green's function representation technique and the corresponding decomposition of the integral kernel. We will concentrate on the solution of 1 + 1 dimensional wave problems with time-independent, non-negative potentials V (x); see the original paper by Leaver [56] and the review articles [61][62][63] for further details. Therefore, we consider the Cauchy problem with t > 0, x ∈ R, and initial data We assume φ 0 and π 0 to be smooth and compactly supported. In this paper we will be interested only in the late time behavior.
The solution of the problem in Eqs. (1) can be repre-sented in terms of an integral kernel k(t, x, y), x, y)π 0 (y)dy, (2) where k(t, x, y) can formally be written as with G(s, x, y) = G(s, y, x) the Green's function, Here f + (s, x) and f − (s, x) are two nontrivial solutions of the mode equation s 2 f − f ′′ + V (x)f = 0 with exponential decay at x → +∞ and x → −∞, respectively, when Re(s) > 0, and the prime denotes the derivative with respect to the variable x. W (s) is their Wronski determinant defined as Note that W (s) does not depend on x, and by definition, it is zero if and only if the two solutions f + (s, ·) and f − (s, ·) are proportional to each other. A priori, the Green's function is a well-defined, analytic function on the right complex plane Re(s) > 0. However, it is interesting to consider its analytic continuation to the left complex plane Re(s) < 0 since this offers the possibility to "close" the integration contour in Eq. (3), as illustrated in Fig. 1. As a consequence, the kernel typically splits into a sum of three different contributions, k(t, x, y) = k modes (t, x, y) + k tail (t, x, y) + k hfa (t, x, y).
(6) Here, k modes is a sum over the residua of the poles of the analytic continuation of the Green's function, and each term in the sum is proportional to e snt . These terms describe the resonances; they oscillate with frequency ω n = Im(s n ) and their amplitude decays exponentially in time since Re(s n ) < 0. Next, k tail is the contribution that comes from the integration around the branch cuts, and it is usually associated with the tail decay, giving rise to a polynomial decay of the form t −p for some p > 0. Finally, k hfa is the contribution from the high-frequency arc (see the curve µ R in Fig. 1). This contribution is highly dependent on the initial data; however, it vanishes for large enough t.
Of course, the structure of the analytic continuation of the Green's function (including the location of the poles and the branch cuts) depends on the potential V (x) and the problem at hand. In the next section, we analyze this structure for the case of a simple toy model, and we show that all the poles are located close to the resonant band i(−µ, µ) and give rise to quasi-bound states. Based on this analysis we expect the structure of the analytic continuation to be very similar in the case of a massive scalar field propagating on a Schwarzschild black hole.
In this work we are mainly interested in the contribution from the poles, that is, k modes , since in our context this is the part corresponding to the quasi-bound states. The poles, s n , of the Green's function are determined by the zeros of the Wronskian W (s); notice that at such zeros the functions f + (s n , ·) and f − (s n , ·) are linearly dependent. Assuming that W (s) has only simple zeros at s = s n , that is, W (s) = α n (s − s n ) + O(s − s n ) 2 , with the mode contribution to the solution from the n'th mode s n is given by with amplitude When the potential V (x) in Eq. (1a) is real, the poles come in complex conjugate pairs, s −n = s * n , and since W (s) and f ± (s, x) are analytic functions in s, the sum φ modes,n (t, x) + φ modes,−n (t, x) is real.
For a generalization of these expressions to more general foliations by spacelike hypersurfaces, see Appendix A.

III. A TOY MODEL EXAMPLE
Before applying the Green's function technique just described to the scalar field problem on a Schwarzschild black hole we consider the following simple toy model potential [33]: consisting of a delta-barrier with positive amplitude A at x = 0 and a step function at x = a > 0, where V (x) jumps from zero to its asymptotic value µ 2 . (Do not confuse the amplitude A in the potential, with the amplitudes A n of the quasi-bound modes.) As it turns out, this model captures all the rough qualitative features of the scalar field problem. In order to construct the Green's function for the potential in Eq. (9), we consider first the region x < 0, where V (x) = 0, such that Notice that f − (s, ·) decays exponentially at x → −∞ when Re(s) > 0, as required. The coefficients a + (s) and b + (s) are unknown to this point and will be determined by the matching conditions. In the region 0 < x < a, the two solutions have the form with coefficients c ± (s), d ± (s). Finally, in the region x > a, where V (x) = µ 2 , we have with coefficients a − (s) and b − (s); and where we have defined κ(s) := s 2 + µ 2 with the choice for the sign such that Re(κ(s)) > 0 for Re(s) > 0, such that f + (s, ·) has the required asymptotic behavior at x → +∞. Because of the delta-barrier, the matching conditions at x = 0 consist of the continuity of the functions f ± (s, ·), and the jump condition f ′ ± (s, 0 + )−f ′ ± (s, 0 − ) = Af ± (s, 0) for their derivatives. At x = a, both f ± (s, ·) and f ′ ± (s, ·) need to be continuous. These matching conditions yield the following values for the coefficients c ± (s) and d ± (s): This information is already sufficient to compute the mode solutions φ mode (t, x) inside the potential well, 0 < x < a, assuming that the initial data is supported inside this well. For this, we first compute the Wronskian, where we have factored out the non vanishing terms c − (s) and c + (s). From this we conclude that the zeros of the determinant are equal to the zeros of the function With these observations Eqs. (8) yield for 0 < x < a, with the amplitude C n given by and In deriving these equations we have used the fact that, for 0 < x < a, Therefore, in order to compute the mode solutions, one first needs to determine the zeros of the Wronskian, which coincide with the zeros of the function F (s) defined in Eq. (15). Next, one needs to determine the quantities B n from Eq. (18), and from this one can compute the amplitude for arbitrary initial data φ 0 and π 0 from Eq. (17). In the following, we compute the poles and the corresponding amplitude coefficients B n in the limit µ ≪ A.

A. Computation of the poles
Before computing the zeros of the function F , it is important to discuss the analytic continuation of the function κ(s) = s 2 + µ 2 for Re(s) ≤ 0. For what follows, we define which is well-defined for all s ∈ C except those for which the argument of the logarithm is negative or zero. Therefore, κ(s) is defined for all s belonging to the set see Fig. 1. This continuation has the property that Re(κ(s)) > 0 for all s ∈ U , implying that the function f + (s, ·) decays exponentially at x → +∞ for all such s. Note that this is different than what occurs for the quasinormal modes, where for Re(s) < 0 the function f + (s, ·) grows exponentially when x → +∞. It turns out that for our purposes the definition in Eq. (20) provides the correct choice of analytic continuation, since it contains the complex interval i(−µ, µ) (the resonant band), close to which the quasi-bound modes lie. For the following calculations it is convenient to replace s ∈ U with the complex angle ϕ defined by the analytic, one-to-one transformation with Q := {ϕ ∈ C : −π/2 < Re(ϕ) < π/2}. The lines Re(ϕ) = const. describe hyperbolas in the s-plane with focal points ±iµ, and when Re(ϕ) → ±π/2, these hyperbolas approach the branch cuts ±i[µ, ∞). The lines Im(ϕ) = const., in turn, describe semi-ellipsis in the splane with focal points ±iµ. In the limit Im(ϕ) → 0 these lines approach the complex interval i[−µ, µ]; see Fig. 2.
Working with the complex angle ϕ instead of s simplifies the following calculations. For example, κ(s) = µ cos(ϕ) and κ ± s = e ±iϕ have simpler representations, and the function F defined in Eq. (15) has the form where the parameter q is defined as q := aµ. As shown in Ref. [33], the zeros of F can be determined analytically based on a series expansion in λ := 2µ/A. When λ = 0, the zeros of F (0, ϕ) are given by whose solutions ϕ (0) n ∈ [−π/2, π/2] are real and satisfy They give rise to a finite set of purely imaginary frequencies s = iω For small values of λ > 0 we introduce the quantity q n := 2q + 2/ cos(ϕ (0) n ) , and the zeros of F (λ, ϕ) can be written as the following series expansion: which yields n ). This expansion shows that, when λ > 0, the poles, which for λ = 0 lie all on the imaginary interval i(−µ, µ), wander to the left half-plane, with the real part scaling like λ 2 for small λ; see quasi-normal mode case, where there are poles arbitrarily far away from the origin.

B. Computation of the coefficients Bn
In order to compute the coefficients B n defined in Eq. (18), we first notice that Here, we can substitute the expansion (26) for ϕ n (λ), which yields sin ϕ n (λ) = sin ϕ (0) cos ϕ n (λ) = cos ϕ (0) It is instructive to discuss in more detail the limit λ → 0, corresponding to a delta-barrier with infinite amplitude, i.e. a hard wall. In this case one expects normal modes, i.e. purely imaginary poles for the Green's function. Indeed, when λ = 0, we saw that the zeros of the Wronskian are purely imaginary, s n ∈ (−π/2, π/2) are the solutions of Eq. (25). In the limit λ → 0, the amplitude coefficients B n reduce to where we recall that ω n . Taking into ac-count that 2s n /A = iλ sin ϕ n (λ), the expression (17) for the coefficients C n reduces to Therefore, the mode solution in this limit is for 0 < x < a, and as expected, it is a sum over purely oscillating factors.

D. Numerical results
To complement our previous analytic results, we solve numerically the Cauchy problem in Eqs. (1) and (9). First, we define a new set of variables, π = ∂ t φ and D = ∂ x φ, to convert Eq. (1a) into a first-order system of linear differential equations. Then, we apply the method of lines with a third order iterated Crank-Nicholson integrator, and use second order finite differences for the spatial discretization. We work on the spatial domain [0, x max ], with x max larger than a. The boundary condition imposed at the origin, where the labels refer to the left and right of x = 0. This condition comes from the behavior of the function dictated by the delta-potential. To see this, notice that due to the continuity at x = 0, φ must satisfy that φ R = φ L , and ∂ x φ R = ∂ x φ L + Aφ L . If we impose that the flux is purely outgoing at x < 0, we have an extra condition, namely ∂ t φ L = ∂ x φ L . After eliminating ∂ t φ L from these equations, we arrive to the desired condition. At the right boundary x = x max we use an outgoing Sommerfeld-like condition.
We find a very good agreement between the numerical solution and the function given by the sum of modes defined in Eq. (16). In Fig. 4, we present the numerical solution and the function obtained from the sum over the contribution of the modes. This example corresponds to a Gaussian initial data of the form φ 0 (x) =φ 0 e −(x−x0) 2 /σ 2 and π 0 (x) = D 0 (x). After an initial transient period, the contribution from the quasi-bound modes dominate and the solution can be accurately represented by the sum of these modes.  16)) and the numerical solution. The parameters for the potential used in this example are: A = 6, µ = 0.3, a = 50, and the parameters for the initial data are:φ0 = 0.1, x0 = 15, σ = 1. The relative numerical error is estimated to be less than 1%.
The fact that at large times the solution can be accurately described solely in terms of the quasi-bound modes is remarkable, since as discussed in Sec. II the solution is really a sum over three contributions. In order to get some understanding about this fact, we analyze in Appendix B the contribution from the high-frequency arc and show that it vanishes identically for times larger than 2a, corresponding to two crossing times of the potential well. Interestingly, it seems that the tail contribution is also negligible, at least for the time scales of our simulations. We will find a similar behavior in the next two sections when discussing the Klein-Gordon equation on a Schwarzschild black hole.
We also calculate the energy of the field inside the well In Fig. 5 we show the evolution of the energy normalized by the energy of the initial configuration. At late times, when the modes contribution is dominant, the remaining energy for this toy model is about 0.6% of its initial value.
FIG. 5: After some transient period in which the field bounces between the potential barriers, the energy drops from its initial value and settles down around a fixed value. This late contribution is given by the sum of the modes, Eq. (16).

IV. QUASI-BOUND STATES AND SPECTROSCOPY FOR THE KLEIN-GORDON EQUATION
The evolution of a scalar field with fixed angular momentum number, ℓ, propagating on a Schwarzschild background spacetime is described by the Cauchy problem in Eqs. (1), where the potential is given now by Here x := r * = r + 2M log(r/2M − 1) is the tortoise radial coordinate, with r > 2M ; please see Ref. [33] for further details concerning the notation in this section.
The physical field Φ is related to the field φ satisfying Eq. (1a) by with Y ℓ a spherical harmonic function with associated angular momentum number ℓ.
Assuming again that all zeros of the Wronskian are simple, we apply the formulae in Eqs. (8) to the Cauchy problem with potential given in Eq. (34) . The difficulty now is that the functions f ± (s, x) and their Wronskian cannot be computed in an explicit simple form, as it was possible in the previous section. As explained in more detail below, the functions f ± (s, x) belong to the confluent Heun class, and it is not known to the authors whether or not there are simple expressions for the Wronskian of such functions. Approximations which are valid in the small M µ limit have been worked out in Ref. [37]. However, in this work we will not make use of these approximation; rather, we will construct the functions f ± (s, x) numerically via a matching procedure.

A. Exact mode solutions
For practical reasons it is convenient to define the dimensionless variables Introducing the ansatz φ(t, r * ) = e st f (z) into Eqs. (1a) and (34), we obtain in terms of these new variables Here we have defined Ω(s) := 2ε 1 + where the function κ(s) = s 2 + µ 2 and its branch points are given in Eq. (20). Notice that accordingly, Ω(s) has positive real part for all s in the domain U defined in Eq. (21).
Eq. (36) has an irregular singular point at z = ∞, and two regular singular points at z = 0 and z = −1, respectively. As z → ∞, the solutions behave as e ±Ω(s)z ; here we are interested in the one with the minus sign, corresponding to an exponential decay. At z = 0 the characteristic exponents of the regular singular point are ±2M s, meaning the solution behaves as z ±2Ms close to z = 0; the plus sign corresponds to the solution we want to consider which is purely outgoing (that is, moving towards the horizon). At z = −1 there is a degenerated characteristic exponent equal to one.
An exact solution to the Eq. (36) satisfying the required conditions at z = 0 is given by f − (s, z) = (z + 1)z 2Ms e −Ω(s)z HeunC(α, β, γ, δ, η, −z), (38) with HeunC the confluent Heun function as defined in MAPLE, the parameters being Since by definition HeunC is regular in the vicinity of z = 0, and equal to one at z = 0, expression (38) gives a closed-form representation for the left solution of our problem.
The quasi-bound states correspond to those frequencies for which the left solution, f − (s, z), decays exponentially when z → +∞. In view of the closed-form expression, Eq. (38), this means that HeunC should not grow faster than polynomially when z → ∞. Unfortunately, the authors are not aware of known asymptotic expansions for the confluent Heun function HeunC when z → ∞.

B. Hydrogen limit
For our applications to dark matter halos discussed below, we are interested in small values of ε = M µ > 0. In this case, the potential V (r) develops a well whose minimum is located at r min ≃ ℓ(ℓ + 1)M/ε 2 . Therefore, focusing on effects that occur near or to the right of the minimum of the potential well, it is appropriate to replace z with the rescaled variable ζ := ε 2 z. Accordingly, we can replace Ω(s) with k(s) := ε −2 Ω(s). When written in terms of ζ and k(s), the mode equation (36) when taking the pointwise limit ε → 0 and assuming that k(s) → k (0) (s) converges to a finite value in this limit. This is the eigenvalue problem encountered in the discussion of the hydrogen atom. Therefore, in the limit ε → 0, the frequencies are determined by the Balmer spectrum for k (0) (s): with the corresponding eigenfunctions Here (a) 0 := 1, and (a) j := a(a+ 1)···(a+ j − 1). In view of Eq. (37), and taking k(s) = ε −2 Ω(s), Eq. (41) implies the following asymptotic expressions for s (compare with Eq. (49) in Ref. [50]): corresponding to a purely imaginary spectrum. This is to be expected, since our considerations here completely neglect effects near the horizon, and as a consequence of our rescaling, the potential is replaced by the effective potential of the hydrogen problem, which diverges at ζ → 0. Therefore, the effect of our zeroth order approximation is to replace the potential barrier by an infinite barrier. Consequently, one finds normal modes instead of quasibound states. This is similar to what occurred in the limit λ → 0 of the toy model problem discussed in the previous section.
As it turns out the zeroth order expression for the frequencies given in Eq. (43) yields already a rather good approximation for the imaginary part of s n , as long as M µ ≪ 1. When ε is nonzero but small, one still expects the solution f to be described by Eq. (40) for large enough z. The solution of this equation which decays as z → ∞ is given by [37] where ξ := 2Ω(s)z, ν = 2/k(s), and U (a, b, ξ) denotes the confluent hypergeometric function as defined in Ref. [64].

C. Shooting to a matching point
Based on the exact expression for f − given in Eq. (38), and the approximate expression for f + given in Eq. (44), we determine by a shooting algorithm the complex frequencies, s n , corresponding to the quasi-bound modes. Starting from a given point r in the interval 2M < r < 4M , and using the known convergent series expansion of HeunC(α, β, γ, δ, η, −z) at z = 0 (with convergence radius one, see Ref. [65]), this algorithm first numerically integrates the mode equation (36) in order to extend the function f − (s, z) to some intermediate point z 1 . (Typically, we choose z 1 = ℓ(ℓ + 1)M/ε 2 to be located close to the minimum of the potential well.) For the numerical integration we use a fourth-order Runge-Kutta scheme. Then, the solution f + (s, z) is extended to z 1 in a similar fashion, starting from a large value of z where the approximation in Eq. (44) can be used. Then, the Wronski determinant, W (s), defined in Eq. (5) is computed at z 1 . The last ingredient of the algorithm consists of a Newton iteration scheme whose purpose is to find the roots of the Wronskian W (s).
In order to illustrate the effectiveness of the method we use similar initial data as the one used in Refs. [32,33].
The radial part of the initial data is given by with the normalization N = [2/(R 1 − R 2 )] 8 . The free parameters R 1 and R 2 allow to set different "locations" and "sizes" for the scalar field distributions at t = 0. The initial value of the time derivative is chosen asφ| t=0 = π 0 (r) = 0.
Next, we show the results for two evolutions corresponding to M µ = 0.30 and M µ = 0.20, respectively. In both cases ℓ = 1, and the initial data is supported inside the interval [R 1 , R 2 ] = [4M, 8M ], see Eq. (45). The first 9 quasi-bound frequencies s n are shown in Table I.
After determining the quasi-bound frequencies, we compute the amplitude each of them is excited from the given initial data. This is performed based on Eq. (A5) of Appendix A, assuming the initial data is specified on an ingoing Eddington-Finkelstein slice. The corresponding amplitudes A n for the first 9 modes are shown in Table II. The number of modes that are excited is infinite in principle, but as can be seen from Table II, the amplitudes decrease rapidly for large n. More specifically, we have found that for the particular initial data chosen, the first five modes are relevant, while the sixth one starts to contribute in a negligible way. Of course, the number of modes requires to accurately describing the modal part of the solution may change depending on the initial data.
The top panels of Fig. 6 show the final output for the excited modes, ψ n (t, r) := Re A n e snt f + (s n , r)/r , where the values of A n and s n are the ones given in Tables I and II. The sum of these modes, ψ(t, r) = n ψ n (t, r) (n = 1, . . . , 9 in our examples) are shown in the lower panels of Fig. 6. In the next section, these solutions will be compared to the numerical evolutions obtained from the same initial data; see Figs. 9 and 10 below.

V. DYNAMICAL EVOLUTION AND COMPARISON
In this section we present some numerical evolutions, with the purpose of extending previous work, as well as comparing with the semianalytical results of Sec. IV. We will corroborate that, after some transient initial period (and for the time scales reached with our numerical evolutions), the scalar field can be accurately described by a superposition of quasi-bound modes alone. Furthermore, the amplitude of each excited mode can be computed using the recipe presented in the previous section; see for instance Figs. 9 and 10, where the agreement between the numerical evolutions and the semianalytical estimations  6.392 × 10 −6 6.90 × 10 −8 3.24 × 10 −8 is remarkable. As mentioned above, only the first few modes contribute significantly for the particular initial data considered in this paper.

A. Dynamical evolutions
Evolutions of quasi-bound modes on a Schwarzschild background were performed in Ref. [32]. In that case, initial conditions were specifically chosen in order to excite mainly individual quasi-bound modes. To our knowledge, long lasting evolutions for more generic initial data, showing the excitation of a combination of quasi-bound states at late times, were first performed in Ref. [33] and later also in [43]. In both cases a generic feature is observed: the excitation of quasi-bound modes, which can be determined (from a Fourier analysis) by a perfect match between the frequencies in the evolution and those of the quasi-bound modes.
We begin this section by extending the study presented in [33] regarding the evolution of arbitrary initial data, considering now a larger number of cases. While the initial configurations evolved in previous works were quite wide (from 100M to 2000M ) as compared to the black hole size, here we consider initial data as thin as 2M .
Starting from the initial configurations described in the previous section, we proceed to solve the dynamic equations to obtain the time evolution. The system is solved numerically with the code used in our previous work, with the addition of (fixed) mesh refinement based on the method presented in [66]. This addition resolved the appearance of noise originated at the right boundary when performing very large evolutions, as reported in [33], by allowing the use of a much larger spatial domain. Aside from this addition, the original code, of which a detailed description is given in [67], was unchanged. To give a brief description, the code uses second order finite differences in space and evolves in time using a method of lines with a third order Runge-Kutta integrator. The use of horizon penetrating coordinates allows one to use free left boundary conditions (inside the horizon), while the right boundary is set far away and all incoming modes are set to zero there.
For the results presented here we set a spatial domain [r min , r max ] with r min = 1M and r max up to about 420, 000M , always satisfying r max > T + 1200M , where T is the total evolution time. This choice ensures that the region [r min , 1200M ], where we analyze the results, is not affected by unphysical signals originated at the right boundary. A total of 8 grids were used, each one doubling in resolution towards the left, with the highest resolution being ∆r = 0.1M . In all cases the finest grid extended from r min to r = 400M , hence covering both the black hole and the initial scalar field distribution. Following Ref. [33], we study the spectrum by evaluating the Fourier transform in time of φ(t, r) at many locations r, and then calculating the average.
The evolutions considered here correspond to initial configurations with scalar field "pulses" with  Fig. 7, was already presented in [33], with the difference that in that work there was some noise to the right of the resonant region (indicated in gray in the figure). In that occasion we claimed that the noise was originated at the right boundary, and that it did not notably affect the spectrum inside the resonant region. Here we were able to eliminate that noise using mesh refinement to push the right boundary far enough, showing that our previous claims were correct.
Given how generically the quasi-bound modes seem to appear in dynamical evolutions, one may wonder how easy it really is to see them manifest as an effective scalar hair (or wig) in general. Put another way, we can ask which conditions allow, and furthermore favor, long lasting scalar field distributions around black holes. To answer this question, we first remind the reader that quasibound modes on a Schwarzschild spacetime cannot exist unless the resonance condition is satisfied [32]. This represents a small range of M µ for typical values of ℓ. For example, for ℓ = 0 the condition is 0 < M µ < 1/4. Second, even when the parameters are in the resonance region, depending on the initial configuration, the quasi-bound modes may get excited with a very small amplitude relative to the initial scalar field distribution. While for very wide initial configurations most of the scalar field remains as a combination of quasi-bound modes, thinner configurations (of order M for instance) lose most of their mass by falling into the black hole, or being radiated towards infinity. This difference between wide and thin initial configurations can be noted by comparing the amplitudes in Fig. 7 (wide initial configuration) with those in Fig. 8 (thin initial configuration). Although the conditions leading to very long lasting scalar fields around black holes are somehow restrictive, we would like to emphasize that they may be satisfied in astrophysically relevant scenarios, as discussed in Ref. [33].

B. Comparison with the semianalytical results
Now we proceed to compare the numerical evolutions with the semianalytical results presented in Sec. IV. For definiteness we consider two evolutions, corresponding to M µ = 0.30 and M µ = 0.20, and in both cases we fix ℓ = 1. The initial data is supported inside the interval [R 1 , R 2 ] = [4M, 8M ], see Eq. (45) above. We compare the results of the evolution both in time and in spatial distribution. In order to compare as a function of time, we fix an observer at a point r = 81.05M , although we have verified that similar results appear for various other values of r. These comparisons are shown in Fig. 9. As also pointed out in Ref. [43], we can see a "beating" due to the combination of different quasi-bound frequencies. We note that, for this particular case, evaluating the contribution from the first nine quasi-bound modes was enough to have a good matching with the dynamical evolution.
Concerning the distribution along the radial coordinate, we can compare the semianalytical result at a given time with a snapshot of the numerical evolution. In particular, we will take the function ψ(t = 104857.6M, r) obtained as the superposition of the first nine quasi-bound modes shown in the lower panels of Fig. 6. The result of such comparison is presented in Fig. 10, where we can see an excellent agreement between the two methods.
To perform a more quantitative comparison it is convenient to first do an error estimation for each case. For a semianalytical solution obtained from the first m modes, Ψ A (t, r) = m n=1 Ψ n (t, r), we estimate the error as the last term in the series. For a numerical evolution with resolution ∆r, Ψ N (t, r) = Ψ ∆r (t, r), we estimate the error as the difference with a solution with resolution 2∆r. In order to do a comparison in time, we take the L 2 norm in r of the difference between the analytical and numerical solutions and their respective errors: We consider that a good agreement up to numerical error is achieved whenever d AN e A +e N . These quantities are shown in Fig. 11 for the same cases shown in the previous figures, that is m = 9 and ∆r = 0.1M in the inner grid. As expected, we see good agreement only at late times, when the semianalitycal solution becomes a good approximation. Note that these error estimations can be considered when deciding where to stop the semianalitical series. For the results presented in this work the error e A is of the order of 1% at the later times shown.
The semianalytic method presented in Sec. IV (and tested against numerical evolutions in Sec. V) offers interesting possibilities. Among them, it allows the evalua-tion of the scalar field at any time t without the need of a full numerical evolution, as long as t is large enough such that the contribution of the high-frequency arc vanishes. 1 This allows us to study efficiently the amount of scalar field that "survives" around the black hole in the form of the superposition of quasi-bound modes that are excited by the initial data. We have referred to these remaining scalar field configurations as Schwarzschild scalar wigs [33]. Furthermore, we can estimate the size of these wigs by estimating the size of each quasi-bound mode. Note that even a narrow initial data, with a width of just 4M , can excite a large number of modes. Most of them, as shown in Fig. 6, have radial extension of hundreds or thousands of M .
In order to obtain a rough estimate on the size of the quasi-bound modes, 2 we can calculate the location of the local minimum, r min , of the effective potential V (r) defined in Eq. (34). Since the actual size is usually much larger [32], this value gives a lower bound for the size of the scalar wig. Fig. 12 shows r min as a function of µ for ℓ = 0 and ℓ = 1, and black hole masses of order M = 10 8 M ⊙ . Larger values of ℓ give larger values of r min . The values shown in Fig. 12 seem to be consistent with galactic halos. For ℓ > 0 we obtain sizes larger than those of the typical visible part of galaxies. For the case with ℓ = 0, although the size seems too small, let us emphasize that r min is usually much smaller than the actual size. Moreover, even if the modes with ℓ = 0 were smaller than the typical size of actual galactic halos, in order to describe a galaxy one would actually have to consider combinations of modes with different values of ℓ for a more realistic scenario, hence yielding larger sizes.
Concerning the amount of scalar field in one of those wigs, we will compute now the total energy, E(t), contained in a slice t = const [32]. For the particular case of the Klein-Gordon equation in ingoing Eddington- Finkelstein coordinates it can be shown that At late times the scalar field can be written as a superposition of quasi-bound states, φ(t, r) = n A n e snt φ n (r), where the mode solutions in Eddington-Finkelstein coordinates φ n (r) are related to the mode functions f − (s n , x) defined in Eq. (38) through the relation φ n (r) = (r/2M − 1) −2Msn f − (s n , z); see Appendix A. Given the properties of the Heun function, it follows that lim r→2M φ n (r) = 1.
Using these identities and integration by parts, it can be shown that the energy in Eq. (48) can be rewritten as Therefore, at late times the total energy of the scalar field can be computed solely from the quasi-bound frequencies and the excitation amplitudes corresponding to a given initial data. For the particular case of the initial data defined in Eq. (45) (with R 1 = 4M and R 2 = 8M ), we have computed the energy E modes (t) given in Eq. (49) for large times, t > 10, 000M , and compared it to the initial energy, E 0 , which was computed from the integral in Eq. (48). The evolution of the normalized energy E modes (t)/E 0 for the same initial data but different values of M µ is shown in Fig. 13.
The following comments on the behavior of E modes (t)/E 0 are interesting.
In contrast to the evolution of the energy shown in Ref. [32], where only one quasi-bound state was evolved and an exponential decay was obtained, for the case of an arbitrary data the behavior of E(t)/E 0 shows transitions from one exponential decaying phase to another one. This is due to the fact that at late times, the solution is a superposition of different quasi-bound states. However, note that the energy E modes (t) is a non-increasing function of time, since the time derivative of the right-hand side of Eq. (49) is equal to −2 times a perfect square. Another interesting property is that for given initial data, the amount of energy of the scalar field is initially smaller but lasts for a larger time when M µ is decreased. The reason for this relies in the fact that when M µ → 0, the real parts of the complex frequencies s n , giving rise to the decay of the solution, become very small, so that the configuration can last for a longer time.
The profile in the energy density ρ nlm (t, r) of each individual mode φ nlm is of the form shown in Figs. 6 and 7, Ref. [32]. However, the energy density of a combination of modes with different frequencies will have nontrivial oscillations, due to its nonlinear dependency on the scalar field. Assuming that galactic halos are described by a combination of these modes, such oscillations might be observable, producing changes in the galactic rotation curves over time.
Let us consider a combination of modes with two frequencies, ω 1 and ω 2 . The energy density will be given by the sum of the energy density of each mode (a quasistatic term), plus a cross term (given the quadratic form of the density as function of the scalar field) with frequency ω r = |ω 1 − ω 2 |.
In general, the maximum possible frequency ω r will then be given by the difference of the higher and lower frequencies. This value will be approximately equal to the upper bound given by the size of the resonance region: is the local maximum (minimum) of the effective potential V (r). In terms of periods, we have the lower bound In order to give some estimations, if we choose M = 10 10 M ⊙ , µ = 10 −22 eV, and ℓ = 0, we have T b = 17 years, thus in this case it seems plausible to detect the oscillations by, for instance, comparing rotation curves observed during a few decades of observation. This relatively small value of T b is obtained for minimum ℓ, and large M and µ, although still within acceptable limits for galactic halos. However, as soon as we choose slightly smaller M or µ, or larger ℓ, the period T b gets extremely larger. For example, M = 10 8 M ⊙ , µ = 10 −23 eV and ℓ = 0 gives T b = 5500 years; M = 10 5 M ⊙ , µ = 10 −22 eV and ℓ = 1 gives T b = 10 11 years. Thus, it remains unclear whether the oscillations might be observed.

VII. CONCLUSIONS
We have considered the late time behavior of localized scalar field configurations surrounding a Schwarzschild black hole, using both numerical evolutions and the Green's function representation technique. At any time in the evolution, the scalar field can be represented as the convolution of the initial data with a suitable integral kernel. The integral kernel can be divided into three different contributions, where each one is a contour integral in the complex frequency plane over the Green's function of the problem. After some transient initial period, and for the time scales reached with our numerical evolutions, we have shown that the scalar field can be accurately described by a superposition of quasi-bound modes alone, which correspond to the residua of the poles of the Green's function.
Given arbitrary initial data for the scalar field on a spacelike hypersurface, the amplitude of each quasibound state in the solution can be obtained by computing a simple, one-dimensional integral. Therefore, once these amplitudes have been calculated, it is possible to predict the late time behavior of the scalar field distribution without performing a numerical evolution.
These conclusions were reached by applying the Green's function technique first to a toy model consisting of a potential formed by a delta and a step function, where the modes can be computed in exact form, and where the asymptotic form of the frequencies and the amplitudes can be obtained analytically. In particular, we showed that for this model problem all the poles of the Green's function are located very close to the resonant band, and thus have very small decay rates. By comparison with a numerical evolution we showed that at late times the solution is accurately described by a superposition of quasi-bound modes. Then, we turned to the case of a massive scalar field on a Schwarzschild black hole background, for which the quasi-bound frequencies, modes and amplitudes were determined semianalytically, and showed that in this case too the agreement with the numerical evolutions at late times is remarkable.
We discussed some of the astrophysical implications of our findings, in particular to the scalar field dark matter models and to the axiverse. We corroborated the results in our previous works, where we found that for appropriate values of the scalar field and the black hole mass the quasi-bound states can be extremely long lived. We also estimated the size of the wigs for supermassive black holes, and speculated about a possible observational feature in the rotational curves of galaxies.
Although the method presented in this paper has been applied only to the case of a massive scalar field living on a Schwarzschild spacetime background, it can be naturally extended to higher order spin fields on other stationary spacetime geometries.
As an example, for ingoing Eddington-Finkelstein coordinates we have with r > 2M , and F (s, x)dx Appendix B: Contributions from the high-frequency arc and the branch cut to the solution of the toy model In Sec. III we provided an analytic derivation of the mode part of the solutions to the toy model problem. In this appendix, we briefly discuss some properties of the remaining parts of the solution, namely the ones coming from the integration over the high-frequency arc µ R and the two curves Γ R around the branch cut (see Fig. 1).
Now that we have control on the terms in the denominator, it is not difficult to estimate each term of the Green's function. For example, for the first term we find e − cos α sinh(R)(τ +ξ−2q) dα.
This concludes the proof of the proposition.