{'fraction_non_alphanumeric': 0.03799034219625611, 'fraction_numerical': 0.012116966734231549, 'mean_word_length': 4.408867392180572, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-gravitating black hole scalar wigs

It has long been known that no static, spherically symmetric, asymptotically flat Klein-Gordon scalar field configuration surrounding a nonrotating black hole can exist in general relativity. In a series of previous papers we proved that, at the effective level, this no-hair theorem can be circumvented by relaxing the staticity assumption: for appropriate model parameters there are quasi-bound scalar field configurations living on a fixed Schwarzschild background which, although not being strictly static, have a larger lifetime than the age of the universe. This situation arises when the mass of the scalar field distribution is much smaller than the black hole mass, and following the analogies with the hair in the literature we dubbed these long-lived field configurations wigs. Here we extend our previous work to include the gravitational backreaction produced by the scalar wigs. We derive new approximate solutions of the spherically symmetric Einstein-Klein-Gordon system which represent self-gravitating scalar wigs surrounding black holes. These configurations interpolate between boson star configurations and Schwarzschild black holes dressed with the long-lived scalar test field distributions discussed in previous papers. Nonlinear numerical evolutions of initial data sets extracted from our approximate solutions support the validity of our approach. Arbitrarily large lifetimes are still possible, although for the parameter space that we analyze in this paper they seem to decay faster than the quasi-bound states. Finally, we speculate about the possibility that these configurations could describe the innermost regions of dark matter halos.

It has long been known that no static, spherically symmetric, asymptotically flat Klein-Gordon scalar field configuration surrounding a nonrotating black hole can exist in general relativity. In a series of previous papers we proved that, at the effective level, this no-hair theorem can be circumvented by relaxing the staticity assumption: for appropriate model parameters there are quasi-bound scalar field configurations living on a fixed Schwarzschild background which, although not being strictly static, have a larger lifetime than the age of the universe. This situation arises when the mass of the scalar field distribution is much smaller than the black hole mass, and following the analogies with the hair in the literature we dubbed these long-lived field configurations wigs.
Here we extend our previous work to include the gravitational backreaction produced by the scalar wigs. We derive new approximate solutions of the spherically symmetric Einstein-Klein-Gordon system which represent self-gravitating scalar wigs surrounding black holes. These configurations interpolate between boson star configurations and Schwarzschild black holes dressed with the longlived scalar test field distributions discussed in previous papers. Nonlinear numerical evolutions of initial data sets extracted from our approximate solutions support the validity of our approach. Arbitrarily large lifetimes are still possible, although for the parameter space that we analyze in this paper they seem to decay faster than the quasi-bound states. Finally, we speculate about the possibility that these configurations could describe the innermost regions of dark matter halos.

I. INTRODUCTION
There is now compelling evidence for the existence of a fundamental scalar field in nature, given the recent discovery of a new particle at the LHC consistent with the Higgs boson in the standard model [1,2]. At a more phenomenological level, however, scalar fields with a deeper substructure have been well-known at least since the dawn of nuclear physics. Nowadays scalar fields are widely used in particle physics, in string theory, and also in cosmology, representing an active area of research in gravitational physics [3].
Along these lines, axion-like particles and other light scalar fields have been considered as possible candidates for dark matter [4][5][6]. At cosmological scales, the highly populated, low-energy excitations of a scalar field could drive the expansion of the universe during the matter dominated era [7][8][9]. Possible mechanisms to excite these long-wavelength modes in a macroscopic way include a regime of thermal equilibrium with the standard model sector giving rise to the appearance of a cosmological Bose-Einstein condensate [10][11][12][13], or a nonthermal vacuum realignment of the axion potential in the early uni-verse [14,15]. Furthermore, it has been shown that the large scale structure that emerges as a consequence of the gravitational instability of the primordial perturbations in this matter component is consistent with observations as long as the mass of the scalar particle is not too light, µ 10 −22 eV [16][17][18][19]. The QCD axion is probably the best motivated and well-known realization of this scenario [20][21][22][23], although there have been many other proposals with a similar spirit [24][25][26][27][28][29][30] motivated to some extent by the absence of any direct or indirect positive signal of the existence of thermal WIMPs in astrophysical observations and/or ground-based detectors.
At the scale of galaxies, the wave properties of an ultralight matter component could in addition alleviate some of the classic discrepancies of the standard cold dark matter model [4][5][6]. On the one hand, the presence of light particles in the universe suppresses the mass power spectrum below a characteristic length scale. For the highly populated zero-mode of a dark matter boson with a mass of the order of µ ∼ 10 −22 eV this scale can be of astrophysical interest, reducing the number of low-mass halo objects from gravitational instability and alleviating in this way the "missing satellites problem" [31]. More-over, using cosmological nonlinear codes for the study of structure formation, large concentrations of ultralight low-energy scalars have been identified to smooth the innermost regions of galactic halos through the formation of solitonic wave-like objects in their centers [32][33][34][35][36]. These boson stars admit a description in terms of a classical field theory [37,38], and more interestingly can alleviate the so called "cusp/core problem" of the standard cosmological scenario [39]. Although possible sources of tension between cosmological and galactic observations in this model have been recently identified in e.g. Refs. [40][41][42], an ultralight scalar dark matter particle is today a hot research topic of modern cosmology.
However, these solitonic coherent macroscopic scalar excitations in the galactic center are likely not alone since most galaxies are expected to host supermassive black holes. It has been known for some time that if a galactic dark matter halo is described in terms of a boson star, then such a configuration cannot live forever [43]. The interaction of an ultralight scalar field with a black hole has been considered in some detail in Refs. [44][45][46][47][48][49][50][51], and extremely long-lived scalar field configurations (referred to as quasi-bound states in the following) have been constructed in the test-field regime of the theory [48][49][50]. Moreover, the distribution of a massive scalar field surrounding a black hole has been shown to give rise to very rich dynamics [52][53][54][55] and to field configurations that fall off so slowly to the point that they may co-exist with the black hole for cosmological times [48,49]. These results encourage further studies of the properties of fundamental fields as a solid candidates to describe the dark matter, providing an alternative to the more usual particlelike description.
In a series of previous papers [48][49][50] we provided a detailed analysis for the propagation of a massive Klein-Gordon field on a fixed, Schwarzschild background. Among other results, we showed that when the Compton wavelength of the field is much larger than the Schwarzschild radius of the black hole, configurations which persist in the vicinity of the horizon for cosmological times do exist [48], that such field distributions are generic [49], and that their late time behavior from arbitrary initial configurations can be described using semianalytic calculations without the need of performing a Cauchy evolution [50].
Clearly, there are several natural extensions of this work towards a more realistic model. One pressing goal is the inclusion of the gravitational backreaction of the scalar field that has been neglected in our previous calculations, implying that our results so far are only applicable to distributions whose mass is small enough such that the gravitational field is dominated by the one generated by the black hole. However, this is clearly not the case for dark matter halos surrounding the black hole at the galactic centers. A further goal is to take into account the angular momentum of the black hole, although the effects of the rotation are probably only important for an accurate description of the scalar field distribution in the region close to the event horizon. Some work in these directions has already been carried out by different groups.
Numerical calculations of the full nonlinear Einstein-Klein-Gordon system of equations for the evolution of massive and massless scalar fields around non-rotating and rotating black holes have been performed in Ref. [53]. The authors report that some fraction of the initial scalar field is accreted by the black hole, leading to an increase in its mass, while another fraction remains in the vicinity of the event horizon with a distribution resembling that of quasi-bound states in the test-field approximation. In Ref. [56], long term evolutions of spherically symmetric scalar field configurations around black holes were also studied by numerical means. Using Gaussian packets as initial data for the scalar field, the evolutions yield configurations describing nonrotating black holes surrounded by scalar clouds. It was found that, for initial configurations consisting of a black hole with a rich scalar environment, after a highly dynamical stage the evolution settles down to states which are consistent with the quasi-bound states. In both works [53,56] it has been found that the scalar field oscillates with a combination of well-defined frequencies, and exhibits a characteristic beating pattern which is due to the combination of two frequencies lying close to each other, as predicted by previous calculations in the test-field approximation [57]. Further numerical work [58,59] has revealed the existence of stationary solutions of the full nonlinear Einstein-Klein-Gordon equations branching off the Kerr solution (dubbed Kerr black holes with scalar hair), provided the rotational parameter of the underlying Kerr black hole assumes some special values. This shows that infinitely long-lived configurations of scalar fields may survive in the vicinity of certain spinning black holes.
The present manuscript provides new results regarding the description of self-gravitating scalar field configurations in the presence of nonrotating black holes, and it extends our previous work in Refs. [48][49][50] to the regime in which the gravitational backreaction of the scalar field becomes relevant. In particular, we present a new method for constructing initial data sets which are devised to represent a time slice of a scalar field distribution surrounding a black hole and surviving the presence of the horizon for extremely long periods of time. Similar states have been previously identified in numerical simulations, but contrary to our construction, these simulations were performed starting from either rather generic initial configurations or from initial data motivated by results from the test-field limit, where large amounts of the field content are lost in the initial stages of the evolution. Furthermore, we also present for the first time an approximate semi-analytic description of these self-gravitating objects, dubbed black hole scalar wigs in this paper, that can be used to explore their characteristic frequencies of oscillation and decaying time scales without the need of performing a time evolution.
In order to achieve this, we introduce an approximation method for solving the Einstein-Klein-Gordon system in the presence of black holes. This technique consists in proposing a coherent state ansatz similar to that of a static boson star, with the important difference that the standard regularity conditions at the origin are replaced by an appropriate "horizon boundary condition", where purely outgoing boundary conditions are imposed for the scalar field. The configurations constructed in this way do not yield exact solutions of the Einstein-Klein-Gordon equations; nevertheless, and as we show, they provide a model for black hole scalar wigs which is accurate for large time spans. Additionally, our configurations can be used to construct exact initial data sets. Based on a numerical evolution of the full Einstein-Klein-Gordon system starting with such sets we prove the viability of our approximation method. The remaining part of this work is organized as follows. In Sec. II we present the Einstein-Klein-Gordon system as well as a simple method for constructing spherically symmetric initial data sets for this problem. Next, in Sec. III we review well-known solutions describing static boson stars or quasi-bound states living on a Schwarzschild background and mention some of their most important properties that are relevant for the investigation that follows. In Sec. IV we present the main result of this paper, namely our approximation method for computing long-lived, self-gravitating scalar field configurations surrounding black holes. This method leads to a nonlinear radial eigenvalue problem which is solved numerically. By construction, these new configurations interpolate between the static boson stars and the quasi-bound states on the fixed Schwarzschild background, and the verification of these limits from the numerical data is also shown in Sec. IV. Next, in Sec. V we perform numerical evolutions of the spherically symmetric Einstein-Klein-Gordon equations without further approximations, starting with initial data corresponding to a fixed time slice of the solutions constructed in Sec. IV. We find that the solutions computed from the numerical time evolution slowly decay in time and oscillate at rates which are in excellent agreement with the approximate solutions, validating our approach. Finally, conclusions and a discussion about the lifetime of some of the black hole scalar wigs models discussed in this article are given in Sec. VI.
Throughout this paper we work with Planck units, in which the gravitational constant, the speed of light and Planck's constant are set to one, such that all quantities are dimensionless. Additionally, we use the signature convention (−, +, +, +) for the spacetime metric.

II. THE SPHERICAL EINSTEIN-KLEIN-GORDON PROBLEM
In this section, we briefly review the equations of motion describing a spherically symmetric, self-gravitating scalar field configuration and introduce the choice of variables and gauge conditions we find useful for our inves-tigation. The main result of this section is the presentation of the Einstein-Klein-Gordon system in the form of Eqs. (5) and (9). As a by-product, we also obtain a simple method of constructing spherically symmetric initial data sets.
Without loss of generality we can parametrize the spherically symmetric spacetime line-element in terms of Arnowitt-Deser-Misner (ADM) variables, Here α denotes the lapse function, β the radial component of the shift vector, γ 2 the radial-radial component of the three-metric, and dΩ 2 = dϑ 2 +sin 2 ϑ dϕ 2 the standard line-element on the unit two-sphere. Due to the spherical symmetry and choice of coordinates, these quantities depend only on time t and the areal radius r, but not on the angular variables ϑ and ϕ. 1 Although the parametrization of the metric tensor in terms of the variables (α, β, γ) is useful for many calculations, in this paper we find it more convenient to work with the slightly different set of quantities (α, ν, m).
Here, as before, the function α denotes the lapse, which determines the future-directed unit normal covector n µ = −α∇ µ t to the slices of constant time t = const., the function ν refers to the quantity which is slice-dependent as well, whereas m is the Misner-Sharp mass function [60] which is invariantly defined in terms of the areal radius r and its differential ∇ µ r as follows: The relation between the ADM variables (α, β, γ) and the new set of quantities (α, ν, m) that we use in this paper is given by The functions (α, ν, m) are related to the matter distribution through Einstein's field equations, G µν = κT µν , which imply In the standard ADM formalism the lapse function and the shift vector are usually associated with the gauge degrees of freedom corresponding to the choice of the spacetime coordinates. In our case, the spatial coordinates are fixed by choosing the areal radius r and the spherical coordinates ϑ and ϕ. Therefore, we do not have the freedom of choosing both α and β arbitrarily anymore; however, we can still choose a combination of these two quantities, as we will do later.
Here κ = 8π is the gravitational coupling constant (remember that we work in Planck units), D 0 := n µ ∂ µ = α −1 ∂ t − ν∂ r is the derivative in the normal direction to the hypersurfaces of constant time, and the prime denotes the derivative with respect to r. The source terms ρ, j and S are defined as the following contractions of the stress energy-momentum tensor T µν with the normal vector n µ and the unit radial vector field w µ = γ −1 δ µ r orthogonal to it: ρ := n µ n ν T µν , j := −n µ w ν T µν , S := w µ w ν T µν . (6) So far these expressions are general and are applicable to any matter fields in spherical symmetry. However, in this paper we will restrict our attention to the macroscopic, coherent, solitonic excitations of a canonical complex scalar field Φ(t, r) with an internal U (1) global symmetry and, for simplicity, no selfinteractions. For this matter component the energymomentum tensor can be expressed in the form Here Φ * denotes the complex conjugate of Φ, Π := D 0 Φ is the conjugate momentum associated with the scalar field, χ := γ −1 Φ ′ , and µ is the mass of the corresponding scalar particle. The Klein-Gordon equation (∇ ν ∇ ν − µ 2 )Φ = 0 for the dynamical evolution of the scalar field can be cast into the following form: Herein, K is the trace of the extrinsic curvature associated with the constant time slices, which can be computed from the functions ν, γ and j = −Re(Π * χ), see below for details. As previously discussed in this section, ν is a gauge parameter which can be freely specified. The lapse function α can then be obtained by integrating Eq. (5c), the radial-radial component of the threemetric γ 2 from Eq. (4), and the mass function m can be computed from Eq. (5a), from Eq. (5b), or from a combination of the two. We end this section with the following important remark. By solving the single equation (5b), we may determine a spherically symmetric initial data set satisfying the constraint equations of the Einstein-Klein-Gordon system. To demonstrate this result, we first note that the components of the extrinsic curvature associated with the hypersurfaces of constant time t are where A, B = ϑ, ϕ. Next, using Eqs. (5a) and (5c) we can rewrite the radial-radial component as Now let m(r) and ν(r) be any two smooth functions satisfying Eq. (5b) with γ = (1 − 2m/r + ν 2 ) −1/2 . Then, we claim that the three-metric γ(r) 2 dr 2 +r 2 dΩ and extrinsic curvature whose components are given by the expressions in Eqs. (11), (10b) and (10c) satisfy the Hamiltonian and momentum constraint equations. To check the claim, we first notice that momentum constraint in spherical symmetry yields the single equation which is trivially satisfied by virtue of Eqs. (10c) and (11). Next, the Hamiltonian constraint is a consequence of the momentum constraint and Eq. (5b), and the claim is then proven. Therefore, given any initial scalar field configuration (Φ(r), Π(r)) and a choice for the gauge function ν(r) one can generate an initial data set for the spherically symmetric Einstein-Klein-Gordon system by solving the single equation (5b) for the mass function m(r). Subsequently, this data may be evolved in time by means of Eqs. (5) and (9) and a suitable gauge condition for ν. We will use this procedure to construct appropriate initial data sets representing a moment of time of selfgravitating scalar wigs surrounding black holes in Sec. IV. These data sets will be evolved numerically in Sec. V using an independent method.

III. BRIEF REVIEW OF STATIC BOSON STARS AND OF QUASI-BOUND STATES SURROUNDING BLACK HOLES
In this section, we provide a brief overview of two different well-known applications of the spherically symmetric Einstein-Klein-Gordon system, Eqs. (5) and (9) above. In both of them the scalar field describes a coherent state of the form with s = σ + iω a complex constant and ψ(r) a complexvalued function depending on the areal radius coordinate r only. Here σ ≤ 0 denotes the decay rate of the state, while ω describes the frequency of its oscillatory behavior. Alternatively, one might think of −is = ω − iσ as a complex frequency with positive imaginary part describing the decay rate. As we discuss next, the ansatz in Eq. (13) is appropriate to describe both static boson stars and quasi-bound scalar field configurations surrounding a Schwarzschild black hole in the test field approximation. Although these solutions have been extensively studied in the literature, we review their most important properties in the next two subsections since they will be essential for the motivation of the more general construction for the self-gravitating scalar soliton wigs of Sec. IV.

A. Static boson stars
Static boson stars are described in terms of globally regular, localized solutions to the spherically symmetric Einstein-Klein-Gordon system, where the field describes a coherent excitation of the form in Eq. (13) with s = iω purely imaginary and ψ(r) real-valued. In the gauge ν = 0, the system in Eqs. (5b), (5c) and (9) reduces to Demanding regularity at the origin and asymptotic flatness at spacial infinity, ψ(r → ∞) = 0, one obtains a nonlinear eigenvalue problem for the frequency ω. Here ψ c and α c are two free positive arbitrary constants. Note however that the system in Eqs. (14) is invariant under the rescaling (α, ω) → λ(α, ω) with some positive constant parameter λ, and hence the value of the lapse at the origin α c is not really physical, although one usually chooses this function in such a way that α(r → ∞) = 1. Solving this problem gives rise to a discrete family of frequencies ω n (ψ c ) for each central value of the scalar field ψ c , where n = 1, 2, 3, . . . labels the different possible solutions. More details about these states are given shortly.
Since the integration of the system is done numerically, in practice one can find the eigenfrequencies ω n (ψ c ) by means of a shooting algorithm. To proceed, one integrates Eqs. (14) outwards, starting from the initial conditions in Eq. (15), setting α c = 1 for convenience, and finetunes the value of the frequency ω to match the asymptotically flat solution ψ(r → ∞) ∼ exp[− µ 2 − ω 2 /α 2 r] at a large but finite value of the radial coordinate r max , increasing the value of r max until the shooting parameter converges numerically. We choose the solution as the one which satisfies the outer boundary conditions for some final r max within a given tolerance. The "correct" values of α c and ω can then be restored a posteriori by making use of the rescaling (α, ω) → λ(α, ω) of the Einstein-Klein-Gordon system, with the value of λ computed such that α(r max )γ(r max ) = 1 (as occurs for the Schwarzschild solution) at the outermost point of the spatial numerical grid. We will follow a similar procedure later in Sec. IV when constructing self-gravitating scalar field configurations in the presence of black holes.
Notice that in order for ψ(r) to decay at infinity, one needs |ω| < µ, which puts an upper bound on the frequencies ω n (ψ c ). The solution corresponding to its lowest possible magnitude, n = 1, describes the ground state, for which the function ψ(r) has no zeros. In contrast to this, the function ψ(r) for the excited states with n > 1 have n − 1 nodes. However, only the ground state is stable under small perturbations [61], and for that reason only this mode can be relevant for late time phenomena.
For a given mass of the scalar particle, µ, the field amplitude at the origin can be varied continuously from zero to ψ crit c ∼ 0.26, increasing in this way the mass m(r → ∞) of the self-gravitating object. Configurations with an amplitude larger than ψ crit c are also unstable against small perturbations [62], and thus they are not of interest for our purposes in this paper. A plot showing the typical profile of the function ψ(r) for the ground state is exhibited in Fig. 1. In this case we have fixed µ = 0.1, although this profile can be easily rescaled to an arbitrary value of the scalar field mass using the invariance of the system in Eqs. (14) under the transformation with ψ, α and γ unchanged. In order to make contact with an ultralight axion dark matter component one would choose a = 8.

B. Quasi-bound states on a Schwarzschild background
Like the case discussed in the previous subsection, quasi-bound states living on a Schwarzschild black hole background are obtained from the coherent state ansatz in Eq. (13). However, contrary to the case of boson stars, we neglect now the self-gravity of the scalar field, which is equivalent to setting κ = 0 in Einstein's field equations. Under these assumptions, integration of Eqs. (5) gives (up to a multiplicative constant which can be reabsorbed in the definition of t) with M BH the (constant) mass of the black hole. The gauge choice ν = 0 leads to the Schwarzschild metric in standard Schwarzschild coordinates, which are regular in the exterior region 2m < r < ∞ but singular at the horizon, r = 2m. The alternative choice ν(r) = 2m/r(1 + 2m/r) −1/2 , which we will use later, leads to ingoing Eddington-Finkelstein (iEF) coordinates, which are regular for all 0 < r < ∞ including the future horizon. 2 Introducing the ansatz in Eq. (13) into the Klein-Gordon equation (9), we obtain, in the gauge ν = 0, Note that the differential equation for the wave-function ψ(r) has precisely the same form as the one in the bo-son star case, Eq. (14c), except that now s = σ + iω has in general a non-trivial negative real part σ < 0, and the metric coefficients α 2 (r) and γ 2 (r) correspond to those of the Schwarzschild black hole, see Eq. (18), instead of being determined by the self-gravity of the scalar excitation. Solutions to Eq. (19) with a purely real frequency −is = ω which decay at spatial infinity still exist; however, they are not regular on the future (nor the past) horizon, and in this sense they do not correspond to physical situations [48]. As mentioned previously, the negative real part σ < 0 is necessary for the solution to decay in time and in this way evading the hypothesis of the no-hair theorems.
As we have discussed in detail in Ref. [50], the quasibound states are obtained by demanding purely outgoing boundary conditions on the black hole horizon, and can be expressed in terms of the confluent Heun function, and HeunC is the confluent Heun function as defined in MAPLE. By definition HeunC is regular close to the origin and equal to one at z = 0, so the expression in Eq. (20) provides a closed-form representation of the scalar field in the vicinity of the black hole horizon. (Note however that since Re(s) < 0, the term z 2MBHs in Eq. (20) is not regular at z = 0, so these scalar field configurations are still singular on t = const. hypersurfaces. We will return to this point in a moment.) For certain discrete values of s, the function defined on the right-hand side of Eq. (20) decays to zero as z → ∞, and these give rise to the quasi-bound states.
Detailed studies of quasi-bound states in Schwarzschild and Kerr black holes have been given in several papers using different techniques, see for instance Refs. [64][65][66][67][68][69]. Perhaps the most accurate way to find the quasi-bound states is via the continued fraction method first employed by Leaver [70] in the context of the quasi-normal modes, and later used by Dolan [67] to find the frequencies of the quasi-bound sates, see also [71,72] and references therein. Furthermore, analytic expressions for the discrete frequencies, valid for small values of the dimensionless parameter ε := M BH µ, were obtained long ago by Ternov et al. [65] and Detweiler [66], and for ℓ = 0 they yield where n = 1, 2, 3, . . .. Although these expressions were only derived for the case with a nonvanishing angular momentum number (after decomposing the field in spherical harmonics), a direct comparison of the frequencies with those obtained from the continued fraction method also yields a very good agreement for the case with ℓ = 0 [49]. These states can be interpreted by reformulating Eq. (19) as a time-independent Schrödinger equation for the field rψ, whose effective potential has a well, see for instance Refs. [48,49] for further details. For the case of interest in this paper where ℓ = 0, the potential only develops a well if the product of the black hole and the scalar field masses satisfies the inequality 3.3 × 10 3 in physical units, so this is a necessary condition for the existence of quasi-bound states if no angular momentum is present in the field configuration [48]. The decay is then due to tunneling of the field through the potential barrier towards the horizon. Although they do probably not form a basis for the solution space, we have shown in [50] that the late time behavior of the scalar field can be accurately described by a linear combination of quasi-bound states, where the coefficients can be computed from the initial data alone by taking appropriate integrals.
The profile ψ(r) of the ground state for the case in which M BH = 1 and µ = 0.1 is plotted in Fig. 2. It is instructive to compare this plot with the one showing a typical profile of a static boson star in Fig. 1. In particular, note there are both real and imaginary contributions to the wave function when the black hole is present. There is again a scaling symmetry similar to that in Eq. (17), with ψ unchanged, that makes it possible to scale the quasi-bound sates to an arbitrary value of the scalar field mass. Note that the product M BH µ is invariant under this symmetry transformation, so we can use the solution in Fig. 2 to reach any value of M BH and µ in the combination M BH µ = 0.1.
For a typical scenario involving an ultralight scalar field one obtains, using the same conversion factor a = 8.2 × 10 −51 µ[10 −22 eV]/µ as in the case of the static bosons stars, as has been previously noticed in Refs. [48][49][50]. However, we emphasize that this expression does not take into account the self-gravity of the boson particles that would constitute the dark matter halo. The correct half-lieftime of these objects which takes into account the gravitational backreaction will be presented in Sec. VI. We end this section with the following important remark concerning the energy associated with the quasibound states. The energy of the scalar field contained in a constant t surface is given by (see Appendix A) where the energy density for the coherent state ansatz in Eq. (13) turns out to be Close to the horizon, the term α 2 |ψ ′ | 2 diverges as z 4MBHσ−1 , and since σ < 0, this implies that E(t) is infinite. Therefore, the test scalar field has infinite energy on Schwarzschild t = const. slices [48,67]. In the self-gravitating case of Sec. IV, such an infinite energy will be clearly problematic, since it implies an infinite change in the mass function, see Appendix A.
On the other hand, when measuring the energy of the scalar field contained in a spacelike slice that penetrates the future horizon, one obtains finite expressions. For instance, choosing iEF coordinates, where the time coordinate t is replaced with the ansatz in Eq. (13) takes the form Φ(t, r) = e stEF ψ EF (r), where in contrast to the function ψ(r), ψ EF (r) = e −Ω(s)z HeunC (2Ω(s), 4M BH s, 0, δ, −δ, −z) is now regular at the horizon, z = 0. The expressions for the metric coefficients, as well as the scalar field energy, are now also manifestly regular, where Applying the balance law (A9) to the solution in Eq. (29), one can simplify the expression for the energy to which shows that E EF (t EF ) decays exponentially as t EF → +∞. This decay is due to the slow accretion of the scalar field mode into the black hole. Correspondingly, along the past direction t EF → −∞, the energy E EF (t EF ) diverges, which explains why the energy E(t) measured on constant Schwarzschild time slices diverges. Clearly, since we are neglecting the gravitational backreaction of these objects, the quasi-bound states described in this subsection only make physical sense if their amplitude is small enough (although at the mathematical level the solutions given in Eqs. (20) and (29) can be rescaled by an arbitrary factor since the differential equation (19) is linear and homogeneous). For large amplitudes, the gravitational effects associated with the scalar field need to be treated self-consistently by solving the coupled Einstein-Klein-Gordon equations. This is the main subject of the next section.

IV. SELF-GRAVITATING SCALAR WIGS SURROUNDING BLACK HOLES
After reviewing the main properties concerning static boson starts and quasi-bound states surrounding Schwarzschild black holes, here we present the main new results of this work, namely the extension of the quasibound states to the self-gravitating regime. After presenting our semi-analytic approximation method for constructing these modes, we provide numerical results and analyze the properties of these solutions. As expected, and as we also show explicitly here, these solutions reduce to the quasi-bound states in the limit in which the ratio between the masses of the scalar field configuration and the black hole is small. In the opposite regime, that is, when the black hole mass is small compared to the total mass in the system, we show that these objects approach standard static boson star configurations. The results of this section will be corroborated with numerical evolutions in Sec. V.

A. The approximation method
Our approximation method is based on the same ansatz in Eq. (13) that leads to the static boson stars and the quasi-bound states of Sec.s III A and III B, respectively. Like in the quasi-bound case but unlike in the boson star one, we now assume s to have a negative real part σ < 0, that is, the scalar field decays exponentially in time. However, unlike the quasi-bound case but like in the boson star one, we now consider the full coupled Einstein-Klein-Gordon system in Eqs. (5) and (9).
In order to obtain the relevant equations in this section, we assume that the scalar field Φ and its canonical momentum Π can be described, approximately, by a coherent state of the form Plugging this ansatz into Eqs. (5b), (5c), and (9), and adopting the ν = 0 gauge yields the following system for the Misner-Sharp mass m, lapse α, and radial wave ψ, functions: with γ = (1 − 2m/r) −1/2 , and where we have used the fact that K = −κrγj/2. These are the equations describing our boson cloud configurations with a black hole in their center, analogous to those in Eqs. (14) for the case of an isolated static boson star (note that indeed these equations reduce to those in Eq. (14) when the frequency is real, −is = ω, and ψ is real-valued). Like in the boson star case, at spatial infinity we still impose asymptotic flatness conditions in order to obtain localized objects with finite mass. However, now the regularity conditions at the origin r = 0 described in Eqs. (15) need to be replaced with appropriate horizon boundary conditions which incorporate the presence of the black hole. This will be discussed next; however, before doing so, it is important to emphasize the following point: the remaining Einstein equation (5a) yieldṡ which shows that the metric cannot be strictly static, since σ = Re(s) < 0. This is also evident from Eqs. (35a) and (35b), due to the presence of the time-dependent factors e 2σt in their right-hand sides. Therefore, our method consists in finding solutions of the system in Eq. (35) in which the factors e 2σt is set to one, for rdependent functions m, α, and ψ, and interpret them as approximate solutions of the Einstein-Klein-Gordon system for which the metric coefficients m and α are time-independent while the scalar field has the timedependent, slowly decaying form in Eq. (13). Although this clearly does not yield an exact solution, it should provide a reasonably good approximation, at least for time scales 0 ≤ t ≪ 1/σ. 3 Alternatively, instead of thinking of the solutions of the system in Eq. (35) as representing approximate spacetime solutions of the Einstein-Klein-Gordon system, we may use them as exact initial data for the full coupled field equations by defining the extrinsic curvature as discussed towards the end of Sec. II. In fact, this is what we will do in Sec. V in order to explicitly check the validity of our approximate spacetime solutions.
After these remarks regarding the nature and validity of our approximation, we now focus our attention to the horizon boundary conditions. Close to the black hole, the boson cloud configurations are expected to have approximately the same form as a quasi-bound state, since in this region the gravitational potential is dominated by the black hole and the metric tensor is expected to be approximately the Schwarzschild one, at least locally. Therefore, we assume that sufficiently close to the inner boundary the scalar field is accurately described by a Klein-Gordon field propagating on a fixed Schwarzschild background of positive mass m 0 = M BH , say, and thus the function ψ(r) should have the form given in Eq. (20). However, as discussed in Sec. III B, this expression diverges as r → r 0 := 2M BH , and further the energy of the scalar field also diverges in this limit, yielding an infinite contribution to the mass function. Therefore, it is not possible to impose horizon boundary conditions at the inner boundary r = r 0 , since both ψ and m diverge there.
For this reason, our strategy is to impose horizon boundary conditions at some radius r 1 = r 0 (1 + z 1 ) slightly larger than r 0 , where the value of the function ψ and its first derivative are determined according to the expression in Eq. (20). More precisely, the horizon boundary condition for the system in Eq. (35) are: with a non-vanishing constant A and where z 1 = r 1 /r 0 −1 and the function f is given by These conditions replace those in Eqs. (15) when a black hole is present. As in the case of boson stars, the value of the lapse function α at r = r 1 is unimportant, because the system (35) with t = 0 is still invariant with respect to the rescaling (α, s) → λ(α, s) by a positive constant λ. Furthermore, since this system is also invariant with respect to a global phase shift transformation, ψ → e iϕ ψ, it is sufficient to take the constant A real and positive. For small enough values of A one should (and one does indeed, as shown further below) reproduce the quasi-bound states, since in this case the selfgravity of the scalar field can be neglected. However, in the method described in this section there is no restriction on the magnitude of A, except that the larger its value, the closer to the black hole horizon we should fix the interior boundary conditions, such that the quantity ∆m(r 1 )/M BH remains small. Now we address the problem of the computation of ∆m(r 1 ). Under our assumption wherein the self-gravity of the scalar field may be neglected close to the black hole horizon, we have seen in Sec. III B that while the energy E diverges on constant Schwarzschild time slices, it is finite when computed on spacelike slices which penetrate the future horizon. Therefore, we can first compute the change of mass ∆m due to the scalar field in iEF coordinates (t EF , r) and then change the result to the gauge ν = 0 we are using in this section.
Given the value for ∆m at some point (t EF , r) = (0, r 0 ), say, the value for ∆m at any other point can be obtained using the integral formula (A12) in Appendix A. Noticing that in the test field limit ψ EF (r) = z −2MBHs ψ(r) = Af (z), we obtain ∆m(t EF , r 1 ) = ∆m(0, r 0 ) + 2δ where we have set δ := κA 2 /4, and where ρ−2M BH j/r is given by Eq. (32) in which ψ EF is replaced by ψ EF (r) = f (z). The integral expression on the right-hand side can be further simplified by applying the balance law (A9) to the solution in Eq. (29), which yields the identity Therefore, the integral in Eq. (39) can be eliminated completely and we end up with the simple expression ∆m(t EF , r 1 ) = ∆m(0, r 0 ) + δ |σ| r 2 0 |s| 2 − e 2σtEF r 2 for the change of the mass function at the point (t EF , r 1 > r 0 ) due to the presence of the scalar field. In the following, we set the integration constant ∆m(0, r 0 ) to zero, which implies that at iEF time t EF = 0 the apparent horizon has areal radius r 0 : 2m(t EF = 0, r 0 ) = r 0 . Notice that in this case, the correction ∆m(t EF , r 1 ) is small for t EF ≥ 0 as long as δ/|σ| ≪ 1 and r 1 lies close to r 0 .
Finally, we transform this formula back to a constant Schwarzschild time slice by setting t EF = 2M BH log(z), such that ∆m(r) := ∆m(t EF , r)| tEF=2MBH log r 2M BH −1 in the horizon conditions in Eq. (37). This concludes the discussion of our horizon boundary conditions, and we are now ready to solve the system in Eq. (35).

B. Results from a numerical shooting algorithm
The system of Eqs. (35) with e 2σt = 1 and boundary conditions as expressed in Eqs. (37) can be solved for any value of s. Nevertheless, since we are mainly interested in self-gravitating compact objects, we need to search for those values of s for which the scalar field decreases asymptotically to zero at spatial infinity, such that the mass m and the lapse α converge to a constant value at large radii. Again, this constitutes a nonlinear eigenvalue problem for the complex frequency −is.
In order to proceed, we follow an approach similar to that outlined in Sec. III A for the boson star case, except that now we start the integration from the inner boundary at r = r 1 where the horizon boundary conditions (37) are imposed, and integrate towards an outer boundary at a large but finite value of the radial coordinate, r = r max , using a second order shooting method. For configurations which are localized in space, it follows from Eqs. (35) that the scalar field must decay exponen- tially at spatial infinity, more specifically, it should behave as ψ(r) ∼ exp[− µ 2 + s 2 /α 2 r] for r → ∞. 4 This behavior is analogous to that of a static boson star; however, now the function ψ(r) and the frequency −is are in general complex. Accordingly, the outer boundary conditions in Eqs. (16) are generalized to This translates into two different conditions for the real, and imaginary parts of the wave function, where λ := µ 2 + s 2 /α 2 , and since s is complex, λ is a complex number as well, i.e. λ = λ R + iλ I . For any given value of the field amplitude A, there is an infinite tower of possible states s n (A) satisfying the system (35) together with the boundary conditions in Eqs. (37) and (42). However, in the following we will focus our attention on the ground state only, that is, the state for which the magnitude of ω is minimal. These 4 More precisely, the asymptotic behavior of the scalar field is given by ψ(r) ∼ (r/2M ) −p e −qr , with p = 1 + 2M q − M µ 2 /q, q = µ 2 + s 2 /α 2 ∞ , and M = m∞ the total mass. However, this only adds small terms of the order 1/r to the right-hand sides of Eqs. (42), (43) and (44), which can be neglected provided rmax is chosen large enough.
self-gravitating objects are referred to as black hole scalar wigs in this paper.
The behavior of the different physical quantities in a typical black hole scalar wig is exhibited in Fig. 3, where the real and imaginary parts of the wave function, ψ(r), as well as the Misner-Sharp mass, m(r), and lapse function, α(r), are shown for a combination of the parameters such that M BH = 1, µ = 0.1, and A = 0.0354. Here the inner boundary has been chosen at r 1 = 2.1M BH , close enough to the apparent horizon, such that the final results are not substantially affected by the choice of r 1 . The scaling symmetry of the Einstein-Klein-Gordon system that has been mentioned previously now implies µ → aµ, s → as, r → a −1 r, with ψ, α and γ unchanged. This makes it possible to map given solutions in the parameter space (M BH , µ, A) to new ones, keeping the product M BH µ and field amplitude A fixed. This feature will turn out to be essential in the discussion section, when making contact with real astrophysical objects. Note that the value of A affects the local solution determining the horizon boundary conditions in Eqs. (37). In particular, increasing the field amplitude gives rise to solutions with higher ratios M T /M BH between the total and black hole masses, where M T is defined as the value of the Misner-Sharp mass function at the final numerical integration point r max , In analogy with static boson stars this behavior is expected at least until a critical value of A is reached, beyond which the ratio M T /M BH starts to decrease and the system becomes unstable. We leave a more detailed analysis regarding this conjecture for future work. Tab. I shows the relevant quantities for some self-gravitating As may be appreciated from Fig. 4, the total mass M T of these configurations with fixed parameters M BH = 1 and µ = 0.1 but varying field amplitude A we have constructed can have their mass as large as several times the black hole mass, so the self-gravity of the corresponding objects is important. For these parameter choices, they have larger decay rates than their associated quasibound states with negligible self-gravity, as can be noticed for instance in Tab. I. The presence of the black hole also compresses the scalar field configurations, making the self-gravitating objects more compact than the corresponding quasi-bound state. This can be seen in Fig. 5, where we compare the energy density of a quasibound state (i.e. no backreaction included), with the corresponding self-gravitating object of same amplitude. For the sake of illustration we also show the profile of a static boson star in this figure, which is regular at the origin.  6: The different limits of a black hole scalar wig. If the mass of the black hole tends to zero, static boson stars are recovered. As long as the field amplitude remains small, the quasi-bound states are a good approximation.

C. Boson stars and quasi bound states as a limit of the black hole scalar wigs
As discussed previously, if the wave function ψ(r) and the mode frequency are chosen to be real, then the nonlinear eigenvalue problem of the previous subsection reduces to the Einstein-Klein-Gordon system describing the static boson stars of Sec. III A. If on the contrary we turn off the self-gravitational interaction of the scalar field, κ = 0, then this same system of equations describes the quasi-bound states of Sec. III B. In what follows, we show explicitly how the black hole scalar wigs smoothly transit between boson star configurations and Schwarzschild black holes dressed with long-lived scalar test field distributions in appropriate limits of the model parameters. (See Fig. 6 for a clarifying picture.) In order to do so, we analyze a family of black hole scalar wigs characterized by the same value of the parameters µ and A, but different black hole masses M BH . In particular, we fix µ = 0.1 and A = 5.66 × 10 −3 , and solve the system in Eqs. (35) as explained in Sec. IV B, starting from M BH = 0.06, and increasing the value of the black hole mass up to M BH = 1. Each configuration in this family has well-defined outputs for the wave function ψ(r) and mode frequency −is, and in what follows we discuss the behavior of these two quantities as function of the total to black hole mass ratio M T /M BH .
In Fig. 7, left column, we plot the real and imaginary parts of s = σ + iω as a function of M T /M BH , for configurations in the above mentioned family. As expected, when the black hole dominates the self-gravitational configuration, i.e. M T /M BH → 1, the black hole scalar wig frequencies approach those associated with a quasi-bound state, in this case the one characterized by M BH = 1 and µ = 0.1. In the opposite limit, when the mass of the black hole is negligible compared to the total mass of the self-gravitating object, M T /M BH → ∞, we recover the frequency of the static boson star with µ = 0.1 and A = 5.66 × 10 −3 . Note that this is a nontrivial result, because even if the black hole is so small that it barely contributes to the overall gravitational field in the exterior region, we cannot think of it simply as a small perturbation of the metric, due to the large compactness ratio 2m/r in the vicinity of the apparent horizon.
Similarly, in the right column of Fig. 7 we show the smooth transition along the aforementioned family of scalar wigs from the boson star configurations to the Schwarzschild black hole dressed with the quasi-bound state, now at the level of the wave function.

V. FULL NONLINEAR DYNAMICAL EVOLUTION AND COMPARISON
In order to validate our approximation scheme and to explore its limitations, in this section we perform numerical evolutions of the spherically symmetric Einstein-Klein-Gordon system. More specifically, we use the selfgravitating approximate configurations constructed in the previous section to provide initial data for the metric and the scalar field, which are then evolved numerically and compared to the approximate solutions. We show that the latter remain accurate for large amounts of time, supporting the validity of our approximation method.
In order to carry out the numerical evolutions, we first rewrite the initial data in terms of horizon-penetrating coordinates. More specifically, we change the time coordinate t associated with the gauge ν = 0 to an iEF one, defined as in Eq. (27), keeping the radial coordinate r unchanged and M BH denoting the apparent horizon mass. In these new coordinates the solutions to our approxima-tion method in Eq. (13) take the form Φ(t, r) = e stEF (r/2M BH − 1) −2MBHs ψ(r), (47) with the function ψ(r) constructed as described in Sec. IV B on the interval [r 1 , r max ]. In order to extend this function on the computational domain [0, r max ] we linearly interpolate it between ψ(0) = 0 and the value it assumes at r = r 1 . Although this interpolation is somehow crude, it occurs entirely within the horizon and thus does not affect the solution in the exterior region.
For the evolution of the initial data we adopt the Baumgarte-Shapiro-Shibata-Nakamura formalism of Einstein's equations in spherical symmetry [73,74], in which the three-metric has the form Here e φ is a conformal factor and a and b are metric coefficients. It has been shown that this formulation is particularly suitable for evolving spacetimes with black holes [56,74]. The gauge conditions we use are the 1+log condition for the lapse and a Gamma-driver condition for the shift vector in order to avoid the slice stretching of coordinates. Although the dynamical gauge conditions adopted here differ from the ones described in previous sections, they have become standard in the numerical evolution of black hole spacetimes and thus we adopt them here as well, providing an extra test to validate the findings in Sec. IV. The initial scalar field distribution described by Eq. (47) at t EF = 0 is shown in Fig. 8. In order to obtain the initial metric and curvature coefficients, we solve the Hamiltonian and momentum constraints in spherical symmetry assuming, without loss of generality, that the initial slice is conformally flat, that is a(0, r) = b(0, r) = 1. Under these assumptions the Hamiltonian and momentum constraints reduce to ordinary differential equations for the conformal factor, e φ , and the radial component of the extrinsic curvature. Furthermore, we write the conformal factor in a puncturelike form assuming that u(r) vanishes at spatial infinity. We solve the coupled equations using a fourth order Runge-Kutta integrator. Initially the shift is set to zero and we used a precollapsed lapse of the form α = e −φ . We numerically evolve the Einstein-Klein-Gordon system for the models A 1 and A 2 of Tab. I. The first case constitutes a small deviation from the test field regime, since the scalar field only contributes little to the total mass of the system, whereas in the second case the mass of the scalar configuration is much larger and comparable to that of the black hole. In both cases we found that the solution in Eq. (47) accurately describes the evolution over large (coordinate) time spans; at least for time periods for which we were able to maintain the code stable (t ∼ 1800) before noise coming from the boundaries interferes with the field.
In Fig. 9 we plot the metric coefficients a(t f , r), b(t f , r) and α(t f , r), where t f ∼ 1800 is the final time of evolution. In Fig. 10 we show the evolution of the maximum (minimum) value of a (b). The behavior of these two fields reflects the evolution of the metric. Initially, the metric is conformally flat and then it changes in time until it settles down to a new state which remains almost unchanged.
A Fourier transform of the maximum value of the real part of the scalar field shows that it oscillates according to Eq. (47). The time of evolution is just large enough to find the values of the frequencies ω 1 = 0.09943 and ω 2 = 0.09914 for the two cases A 1 and A 2 , respectively, that are in a good agreement with the values given in Tab. I.
In Fig. 11 we also monitor the mass of the scalar configuration throughout the evolution by integrating over the density given in Eq. (8a) with n µ = α −1 (δ µ t − βδ µ r ) the four-velocity of the Eulerian observers in the 1+log slice. We found the time decay rate via a numerical fit for this   mass of the form e 2σnt . Our results indicate that after an initial stage of complex dynamics, the mass presents an exponential decay, indicating that the field is falling into the black hole. For model A 1 we find a time rate decay σ n for the field of the order −1.5 × 10 −5 , whereas for model A 2 we obtain −2 × 10 −5 , which are roughly the values of the predicted rates given in Tab. I. These results are consistent with the expression for the time evolution of the energy given by Eq. (33) in iEF coordinates. This result is interesting because it shows that once we have constructed the black hole scalar wigs on a regular foliation (the t EF = 0 slice) it remains finite in the other regular foliation determined by the 1+log slicing condition.
The main constraint in getting longer evolutions is the fact that the time step is limited by the Courant condition which requires it to be sufficiently smaller than the spatial grid size, which in turn needs to be small in order to keep the numerical solution stable. Thus, in our current setup, it is not possible to obtain evolutions for very large times with a spatial domain that is large enough such that the solution is not contaminated by boundary effects.

VI. DISCUSSION AND CONCLUSIONS
In this paper we have extended the quasi-bound states to the regime where the gravitational backreaction of the scalar field is relevant. In order to do so, we introduced a new approximation method for semi-analytically solving the Einstein-Klein-Gordon system in presence of black holes. This method consists in proposing a coherent state ansatz for the scalar field which, together with appropriate (inner) horizon boundary conditions and asymptotic flatness at spatial infinity, leads to a nonlinear eigenvalue problem which is solved numerically. The viability of our approximation has been confirmed using numerical evolutions of the full nonlinear equations. Although in this work we have focused on nonrotating systems, it should in principle be possible to generalize it to more sophisticated scenarios including rotation.
Our results corroborate that, as expected from previous numerical dynamical evolutions [53,56], the quasi-  Fig. 7 for the magnitude of the decay rate |σ| as a function of the total to black hole mass ratio, now in a logarithmic scale. This plot makes the powerlaw behavior for MT/MBH 3 of this function evident.
bound states possess analogues in the self-gravitating case, yielding self-consistent, asymptotically flat solutions of the Einstein-Klein-Gordon equations. We called these self-gravitating objects black hole scalar wigs, in analogy with the Schwarzschild scalar wigs described in our previous work [49,50].
These objects are characterized by three independent parameters: the scalar field mass µ (like is the case for any solution of the Klein-Gordon equation), the black hole mass M BH (in analogy to the quasi-bound states), and a field amplitude A (which plays a similar role to the central value of the scalar field in the static boson star case). 5 However, thanks to the scaling symmetry described in Eq. (45), it is possible to eliminate one of the mass scales and characterize black hole scalar wigs solely in terms of two parameters. It is convenient to express these two parameters in terms of the following dimensionless quantities: the product M BH µ between the black hole and the scalar field masses (in Planck units), and the ratio M T /M BH between the total and black hole masses.
As expected, the black hole scalar wigs change considerably as M T /M BH increases, undergoing a transition from a quasi-bound state in the test field limit when M T /M BH = 1, to cases in which the mass of the scalar wig is several times larger than the mass of the black hole and represents a large boson star with a tiny black hole in its center (the black hole being tiny when compared to the size of the whole self-gravitating object). In the context of the ultralight axion dark matter models these latter situations may represent the large core of the halos in presence of central supermassive black holes.
In order to make contact with real astrophysical systems, we conclude with some rough estimates regarding the lifetime of the scalar field configurations constructed in this article. In order to do so, we first note that the half-lifetime of the black hole scalar wigs can be estimated from the results given in Sec. IV, together with the expression in Eq. (23). Next, as a specific example we consider the family of black hole scalar wigs characterized by µ = 0.1 and A = 5.66 × 10 −3 . Fitting the magnitude of the decay rate as a function of the total to black hole mass ratio to a power law, see Fig. 12, one obtains the empirical formula where a 1 = 1.8582 × 10 −5 , a 2 = −1.9489.
Introducing Eq. (50) into the expression (23) yields The black hole mass M BH in Planck units belonging to a particular configuration in the above family can be read off, or extrapolated, from the upper horizontal axis in As an example, considering a case where the scalar field configuration is hundred times higher than the black hole mass, Eq. (51) yields a time of the order of ten million years, corresponding to a black hole of mass 10 9 M ⊙ . Comparing this with the order 10 13 years obtained from Eq. (24) for a test field configuration with the same value for the product of M BH µ ∼ 10 3 , we notice that the selfgravitating configurations may have much smaller lifetimes, although we stress that this conclusion is based on a very specific family of black hole scalar wigs. We leave a more systematic analysis of these states and their applicability to real astrophysical systems for future work. In practice, the self-gravitating scalar wigs we have discussed in this article may be detected through their dynamical influence on the stars, gas clouds, light rays or other compact objects surrounding the supermassive black hole at the center of galaxies. Indeed, the motion of these bodies is affected by all form of mass contained inside their orbital radius, that is, by the mass of any (hypothetical) matter distribution in addition to the black hole mass. As we have shown, for the case of a scalar field, such distribution can be several times larger than the actual mass of the black hole. Near-future observations such as the Event Horizon Telescope project [75] should yield accurate estimates for the size of the black hole shadow in Sagittarius A * and at the center of other galaxies, thus providing an estimate for the actual black hole mass. If there was a difference between this mass and the total mass measured through independent, dynamical observations, this would be a clear indication for the presence of large matter distributions surrounding the black hole, such as the ones discussed in the present work. A further possible signature for the presence of a scalar wig could be detected in the gravitational wave signal emitted by the accreted matter in the form of a monopole black hole ringdown [53]. However, we leave these and other studies for future investigation.

Acknowledgments
This work was supported in part by CONACYT  surfaces and the application of Gauss' theorem leads to the "conserved charge" where Σ denotes one of the Cauchy surfaces and √ g Σ d 3 x the induced volume element on it. (Here, we have assumed that the matter fields fall off sufficiently fast in the asymptotic region for the boundary terms to vanish. If T µν does not fall off sufficiently fast or if the spacetime region over which Eq. (A2) is integrated is enclosed between two spacelike hypersurfaces Σ 1,2 which are not Cauchy surfaces, then the difference between Q X [Σ 2 ] and Q X [Σ 1 ] is given by the boundary flux integrals, see further below for an example.) In particular, if X µ is a future-directed timelike vector field satisfying S = 0 and if T µν satisfies the dominant energy condition, 6 then the conserved charge Q X is guaranteed to be nonnegative, and in this case it may be associated with "mass" or "energy".
A particular scenario for which S = 0 occurs when X µ is a future-directed timelike Killing vector field, in which case ∇ (µ X ν) = 0 and the quantity Q X in Eq. (A3) is just the conserved energy E associated with the timeinvariance symmetry of the spacetime generated by X µ . However, dynamical spacetimes such as the ones analyzed in this work do not admit any timelike Killing vector fields and in this case one does not expect a priori the existence of a timelike vector field X µ satisfying ∇ µ J µ X = 0. Interestingly, spherically symmetric spacetimes provide an exception to this expectation. In this case, even though the spacetime might be dynamical, there exists a natural choice for a vector field X µ which does lead to a conserved current density ∇ µ J µ X = 0. This vector field is known as the Kodama vector field [76], see also [77]. In order to describe it, we assume the spacetime to be oriented and denote by η αβγδ the natural volume form induced by the spacetime metric g µν . Denote by {e µ 2 , e ν 3 } a local oriented orthonormal frame tangent to the invariant two-spheres, and introduce the anti-symmetric tensor fieldη αβ := η αβγδ e γ 2 e δ 3 . Thisη αβ can be identified with the volume form on the two-dimensional manifoldM orthogonal to the invariant two-spheres. Then, the Kodama vector field is defined as K µ :=η µν ∇ ν r. (A4) By construction, the integral curves of K µ leave r invariant and thus K µ is timelike in the region where the 6 The dominant energy condition states that for any futuredirected timelike vector field u µ the associated current density J µ u = −T µ ν u ν must be future-directed timelike or null [63]. gradient of r is spacelike, that is, in the region outside trapped spheres. As a consequence of Einstein's field equations it is not difficult to verify that K µ satisfies such that the Kodama vector field K µ is a Killing field whenever T µν = 0, or more generally, whenever the projection of T µν onM is proportional to the induced metric onM . In terms of the ADM parametrization of the metric tensor discussed in Sec. II, we find which yields the timelike Killing vector ∂/∂t in the Schwarzschild case. Remarkably, even when K µ is not a Killing vector field, it still yields a conserved current density when contracted with T µν . Indeed, using Eq. (A5) one finds and thus one obtains a conserved charge Q K of the form given in Eq. (A3). For asymptotically flat, spherically symmetric spacetimes, this charge coincides with the total (ADM) mass M ADM of the spacetime, as we show next.
Combining Eqs. (A9), (A10) and (A11) we also find the identity m(t 2 , r 2 ) = m(t 1 , r 1 ) + κ 2 t2 t1 r 2 T µν K µ (γνn ν + w ν )αdt which can be used to determine the mass function at an arbitrary point (t 2 , r 2 ) of the spacetime manifold from its value at a fixed reference point (t 1 , r 1 ) and an integral over the stress-energy tensor. In particular, when t 2 = t 1 = t the identity (A12) reduces to m(t, r 2 ) = m(t, r 1 ) + κ 2 r2 r1 r 2 T µν K µ n ν γdr t , (A13) and as long as the dominant energy condition is fulfilled and K µ is timelike it follows that m(r 2 , t) ≥ m(r 1 , t) for r 2 ≥ r 1 , that is, the mass contained inside a spherical shell is positive. For an asymptotically flat spacetime the Misner-Sharp mass function m(t, r) converges to the total mass M ADM when r → ∞ along a hypersurface t = const. which reaches spatial infinity. The identity (A12) is key to our discussion in Sec. IV. Finally, we note that for the particular case where the matter field is described by a canonical complex scalar field Φ with a U (1)-symmetric potential, there is an additional conserved current density which is associated with the internal U (1) symmetry of the field and is given by The corresponding conserved charge Q defined by Eq. (A3) is usually identified with the difference between the number of particles and antiparticles in the configuration. For the particular case of the static boson stars discussed in Sec. III A this quantity is given by It should be interesting to analyze the implications of the conserved quantity Q for the scalar field configurations surrounding a black hole constructed in this article.