interplay between black holes and ultralight dark matter: analytic solutions

Dark matter (DM) can consist of a scalar field so light that DM particles in the galactic halo are best described by classical waves. We investigate how these classical solutions are influenced by the presence of a non-rotating supermassive black hole at the center of the galaxy, using an analytical, albeit approximate, approach. Relying on this analytic control, we examine the consequences of imposing causal boundary conditions at the horizon, which are typically overlooked. First, we examine the scenario where the backreaction of dark matter can be neglected. The scalar field decays like a power law at large distances, thus endowing the black hole with “hair”. We derive solutions for the field profile over a wide range of parameters, including cases with rotating dark matter. As a by-product, we extract the dynamical Love numbers for scalar perturbations. Next, we determine the spectrum of bound states and their behaviour. Finally, we incorporate the self-gravity of the scalar field, with a focus on the situation where dark matter forms a soliton (boson star) at the center of the galaxy. We derive an analytical expression for the soliton at every distance from the center. With a solution that remains applicable even at horizon scales, we can reliably compute the accretion rate of the black hole.


Introduction
Quoting Edward Witten [1]: "One side of what theoretical physicists do is to try to understand the fundamental equations of nature, and the other side is to try to solve the equations in different situations and work out predictions for what will happen." While this paper discusses a new hypothetical constituent of the Universe in the form of an ultralight scalar boson, our focus will be on the latter aspect.The equation describing this boson is actually quite simple: a massive scalar field with no self-interactions, minimally coupled to gravity.However, our specific interest lies in a particular class of solutions that describe a self-gravitating, spherically-symmetric, stationary configuration of the scalar in the presence of a black hole (BH) at the center.Furthermore, we will adopt an analytical, albeit approximate, approach to study such solutions.
The main motivation for considering ultralight scalar fields arises from the existence of dark matter (DM).Within a wide mass range of 10 −22 − 10 −21 eV to 10 eV, if we assume that all of the dark matter originates from the light scalar and take into account the density of dark matter in the galactic halo, the occupation number of the particles becomes so large that they behave as oscillating classical fields.Hence, this scenario is commonly referred to as wave dark matter (WDM) [2][3][4][5][6][7][8][9][10][11][12][13][14].The most compelling example within this category is the QCD axion, followed perhaps by axion-like particles (ALP) predicted by string theory [15][16][17][18].
At the lighter end of the mass spectrum, the scalar is known as fuzzy dark matter (FDM).In particular, FDM was proposed to address certain issues pertaining to the small-scale features observed in simulations conducted using conventional cold dark matter models [7].
In the presence of gravity, a free scalar field can support itself against gravitational collapse by the quantum property that particles cannot be localized beneath distances of the order of their Compton wavelength.The resulting (quasi-)stationary solutions of the field equations of motion are commonly known as boson stars (for a complex scalar field) [19][20][21][22][23], oscillatons (for a real scalar) [24][25][26][27][28], or solitons (the term we will use in the following).There exists an upper limit to the mass of a soliton M sol , as a function of the scalar field mass, in order to prevent collapse into a black hole.Simulations [12,29] indicate that fuzzy dark matter halos typically exhibit a solitonic core surrounded by a halo with a Navarro-Frenk-White (NFW) profile [30].Numerous papers have compared this expectation against observational data [31][32][33][34][35][36].However, our work deviates from this approach and instead addresses a different question: how is the wave dark matter soliton in a galactic halo influenced by the presence of a supermassive black hole (SMBH) at the center of the soliton?For this reason, we will take a broader perspective than fuzzy dark matter and consider the mass µ of the scalar as a free parameter within the WDM range.
An important element, once additional matter in the form of a black hole is added to the scalar condensate, is the boundary condition at the black hole horizon.In much of the literature, which primarily focuses on distances significantly larger than the black hole horizon r BH , its effect is often neglected.While the causal solution should be outgoing (entering the horizon), in principle the time-reversed mode also exists.When solving the equations from infinity, one generically encounters a superposition of both modes at the horizon.Prohibiting the non-causal mode of the scalar field could impact the solution, potentially even at distant locations.Numerical control is challenging due to the presence of the horizon and the highly oscillating nature of the solution.Therefore, as mentioned earlier, we will instead employ an analytic approach.For other investigations on how a black hole modifies the soliton see [34,[37][38][39][40][41][42].Black hole accretion of diffuse scalar is studied instead in [43][44][45].
The structure of the paper is as follows.We begin Section 2 discussing the case where the backreaction of dark matter can be neglected, and the equation for the scalar field reduces to the Klein-Gordon equation in a Schwarzschild background.In such cases there exists more literature, dating back to Starobinskii [46] and briefly summarized in [14], and the exact solution in terms of a special function, known as the confluent Heun function [47], has been derived.However, manipulating this exact solution is challenging.Therefore, a systematic exploration of the resulting hairy black hole solutions was conducted in [44] for different values of the dark matter Compton wavelength ℏ/µ compared to r BH , in the case where dark matter has no angular momentum.Contrary to solitonic solutions, these profiles decay like a power law at large distances, thus they endow the black hole with "hair" that is then matched to the halo.Notice that no-hair theorems do not apply because the scalar field is not static.
We will address this scenario using different techniques that provide better control over the approximate solutions and allow for the inclusion of spinning dark matter.In the large mass (or particle) limit, µr BH ≫ 1, the method of uniform approximations [48], which generalizes the standard WKB approach, will be employed.In the opposite (wave) limit µr BH ≪ 1, where both WKB and uniform approximations fail, the approximate solution will be derived by employing the technique of boundary layer theory.The particle limit solution will be sensitive to the boundary condition at the horizon at all distances.On the other hand, in the wave limit, the boundary conditions become negligible within a few black hole radii.In addition our solution enables us to compute the dynamical response coefficients (Love numbers) for a massive scalar field in the small mass limit.Where they overlap, our results are consistent with [49,50].
The intermediate mass range µr BH ∼ 1, with non-zero angular momentum, poses additional challenges and will be thoroughly discussed in Appendix A. This range is particularly intriguing as it pertains to cases where the Compton wavelength of the light scalar field is comparable to the black hole horizon, resulting in a faster growth rate of the superradiance instability.To investigate this instability, our computations would need to be generalized to a Kerr background.However, even for non-spinning black holes, our analytic formulas can be compared with the results obtained in [51] by numerically evaluating the exact Heun function solution.
Finally, we end Section 2 by discussing solitons (of mass M sol ) at the center of the dark matter halo whose dynamics is dominated by the SMBH, M BH /M sol ≫ 1.These solutions are similar to the previously discussed hairy black hole profiles up to a length scale r ≲ r BH /|γ|, where µγ is the binding energy of a single scalar particle.At greater distances, the solution is exponentially suppressed.A second difference is that the frequency has a tiny imaginary part related to the inverse time of absorption of the soliton, which is meta-stable.We will show that causal boundary conditions at the black hole horizon impose the upper bound µr BH ≪ 1 for any black hole dominated soliton.When this bound is satisfied, causal boundary conditions do not need to be explicitly checked.
More interesting is the case where the soliton dominates the dynamics, at least far from the SMBH.This can happen in the opposite limit, M BH /M sol ≪ 1, and it will be discussed in Section 3. The soliton mass can be related to the halo mass via the numerically extrapolated soliton-halo relation, however we will simply keep the soliton mass as a free parameter since we do not model the halo.Neglecting the halo will not significantly affect the solution, as long as the soliton dominates over the halo density.In this scenario a natural question arises: can the soliton dominate the dynamics at all length scales, even down to the size of the black hole's event horizon?As we will see, our analysis reveals that this is never the case.Regardless of how small M BH is, there always exists a length scale, denoted as r e , where the effect of the black hole becomes parametrically equal to that of the soliton.Notably, this length scale satisfies the condition r BH ≪ r e ≪ r sol , where r sol corresponds to the characteristic size of the soliton.For r ≫ r e the self-gravity dominated solution was studied both numerically [12,29,41,42] and analytically [52], while for r ≪ r e it reduces to the one discussed by Hui et al. [44] for dynamics dominated by the SMBH.The two solutions will be matched at r ∼ r e by imposing continuity.We will discover that soliton stability under gravitational collapse imposes the upper bound µr BH ≪ MBH M sol (≪ 1) for any self-gravity dominated soliton.We will conclude that also these kinds of solitons are only possible in the (very) small mass regime µr BH ≪ 1, when causal boundary conditions are unimportant.
Conventions: We will set ℏ = c = 1 throughout the paper, together with r BH = 2GM BH = 1 within Section 2. However, we will restore this value in a few key expressions and in section 3 for the sake of clarity.

Black hole domination
Our starting point is a complex scalar field of mass µ in a black hole background.Minimal coupling leads to the equation of motion where the background metric is strictly Schwarzschild and gravitational backreaction is assumed to be negligible, postponing the discussion of a self-gravitating scalar.
Exact solutions to this equation are known in terms of the confluent Heun function.However, these solutions are not transparent or easily understandable 1 .The confluent Heun function is a special function that can be quite complex, making it difficult to gain insight into the behavior of the solutions just by examining the function itself.To give an example, the exact solution does not provide an explicit relation between the density of dark matter at large distance and the density near the horizon.Progress in this direction was achieved by Bonelli et al. in [53], where they express these connection formulas in terms of a series.While their expressions are in principle exact, our goal is to find instead an approximate solution for the scalar background, valid at all distances, which, on the other hand, is expressed in terms of significantly simpler formulas.One of the results of our analysis is the identification of the different length scales that affect the field profile as a function of distance.These length scales determine where the field's behaviour undergoes significant changes.However, these length scales are not readily apparent when considering only the Heun function, emphasizing the need for alternative approaches with simpler expressions.Finally, we extend the analysis of Hui et al. [44] by discussing also the case of spinning dark matter.
The first step is to specify the time and angular dependence of the field.We look for solutions 1 Besides, software like Mathematica has some difficulties in handling and plotting this function.
with definite frequency ω and angular momentum l, m.Restricting to one frequency is particularly natural because the equation is linear.Barring bound states, the field must have energy at least equal to its own mass, and as a first approximation we will neglect excited states thus setting ω = µ.The implication is that the energy density of the field ∼ |∂ t ϕ| 2 is independent of time and so is the backreaction on the metric (which for now we take to be negligible).
A real scalar field would obey the same differential equation but its time dependence would be ∼ cos(ωt).The energy density ρ ≃ 1 2 φ2 + µ 2 ϕ 2 would not vary with time, but the other terms in the energy-momentum tensor would, leading to a small mixing with higher harmonics and a slow decay of the field due to gravitational radiation.

The usual separation of variables ansatz
brings 2.1 to the form (putting r BH = 1) 3) The form of the derivative operator that acts on ϕ suggests the change of variable because the derivative term becomes d 2 ϕ dr 2 .The equation takes the form of a 1d Schrödinger equation where r goes from −∞ being the horizon to r = 0 − being large physical distances.r will only play a formal role in what follows: when talking about distances, we will always mean r.This definition will prove convenient for our purposes because the derivative term of 2.3 will not change once angular momentum is included.We stress that r is not the tortoise coordinate usually employed to analyze this kind of problems.From our perspective, this change of variables is more natural because it is the standard trick that one uses to eliminate the first order derivative in a second order differential equation.This step will allow us to employ WKB methods to look for a solution of 2.3.
Most of the literature assumes that, if we are only interested in r ≫ r BH , we can replace (r − 1) with r in equation 2.3 and the horizon can be entirely neglected.However, as pointed out by Hui et al. [44] and Baumann et al. [54], the role of the horizon scale is two-fold.Besides trivially entering equation 2.3, r BH is also a boundary where we have to impose boundary conditions on ϕ.
To be more explicit, we will see that the two linearly independent solutions of 2.3 are waves decaying at infinity, one infalling and one outgoing.One could imagine setting up a scattering experiment by choosing the amplitude of the infalling wave at infinity, the reflected outgoing wave being completely fixed.What picks the amplitude of the outgoing wave?The answer is the causal boundary condition we impose at the horizon.
Without this boundary condition, the solution will generically be inaccurate at all distances (if one includes other interaction, such as self-gravitation of dark matter, the solution will be affected only if the black hole is dominating the dynamics).This absence of decoupling between black hole horizons and large scales is not unheard of: it is reminiscent of the black hole Love numbers which, although defined in terms of the behaviour of the field at large distances, are tuned to zero by the presence of the horizon.We will actually compute such tidal response coefficients in equation 2.36.Our goal in the second part of the paper will be to check the effects of this boundary condition and develop approximations for the field profile at different distances.

Large µr BH regime
When attempting to find a solution for a challenging second-order differential equation, a very natural approach is to try the WKB approximation.Despite its high accuracy, the WKB approximation is known to be ineffective near turning points, which are locations where the particle would classically come to a halt.To address this, we begin by determining the range of the parameter space (µ, l) where the WKB approximation remains valid throughout.This will translate into a lower bound on µ.
Looking at equation 2.3, the turning points are given by V (r) = 0. We then begin our discussion by studying the sign of V as a function of r.First, are both positive.However, V (r) can develop two zeros at r > 1 and turn negative between them.Solving for V (r) = 0 we get Thus there are turning points for small µ and non-zero l.To be more precise, only the mass range is free of turning points.The location of these points is well approximated by in the small mass regime.We now assume µ ≫ µ c and apply the WKB approximation, leaving the small µr BH case for later.

Uniform approximations
Assuming µ ≫ µ c (l) we could write down the WKB solution of equation 2.3 in the absence of turning points, but this would turn out to be inaccurate for l = 0. To have proper control of the approximated solution, we will rather derive the WKB approximation and the condition for its validity.Another motivation for taking this route is because it allows us to introduce the framework of uniform approximations (which vastly generalizes WKB) in a simple setting.We will use this more general framework in appendix A, where we will deal with two turning points of the potential close to each other.A comprehensive reference for uniform approximation techniques is [48].
The key insight of uniform approximations is to substitute 2.3 with a 'similar' differential equation for a new function f (x), defined over a new domain x ∈ (1, ∞).We want the new differential equation to be explicitly solvable and to have the same singularity structure (i.e.zeros and poles inside the domain of interest) as the original one 2.3.We will call this the 'similarity conditions'.Under these conditions, we will be able to write down an approximate solution of 2.3 in terms of f (x) by suitably "stretching" x into r.
We choose f (x) to obey with Γ(x) ≡ 1. 2 To simplify the derivative term, we define (similarly to 2.4) Since the derivative term has the same zeros in both 2.3 and 2.9, and both V (r) and Γ are always positive, we can expect the solutions of 2.9 to be a slight deformation of those of 2.3.
We formalize this idea by introducing an unknown function x(r), which will quantify this stretching ϕ(r) ≡ dx dr Upon substitution in 2.3 and making use of 2.9, we obtain an equation for x(r) (2.12) where f has dropped out and the second term is a schwarzian derivative we wish to neglect.This term blows up whenever the singularities of the two differential equations don't match, justifying the similarity condition requirement.
We are thus led to a simple differential equation for x(r) together with a condition for the validity of the approximation For l = 0 we can easily solve this, getting x(r(r ) for large r and µ (ln(r − 1) + 2 − 2 ln 2) + O(r − 1) when r ≃ 1.
Given the solutions e ±ix of 2.9 and 2.11, we obtain the sought field profile where the ↔ symbol indicates that we matched the solution close to the horizon with the one at large distance.We emphasize that this matching is straightforward because the approximation is valid in the whole domain.
Subleading terms in the r → ∞ expansion become comparable with the leading contribution at r = 1.Causality at the horizon then selects the purely infalling wave (r − 1) −iµ , thus enforcing the behaviour e −2iµ √ r r 3/4   at large distances3 .This solves the connection problem completely.
Let's draw some conclusions.We observe that, for all l (up to subleading WKB orders), the amplitude of the profile is completely controlled by the overall factor (dx/dr) − 1 2 = V − 1 4 which has a very simple analytic expression.All complications reside in the less important phase of the wave.
We now compute x for l > 0. The integration of 2.13 gives which is hard to evaluate analytically.While a numerical solution is straightforward, let us observe that when l = 0 we are back to the previously studied case, up to a correction which is small uniformly when r ≥ 1. Indeed max l(l+1)r(r−1) 2 , and quickly becomes ≪ 1 in the large µ regime.Thus, upon expansion, r l(l + 1) Higher orders are as just easy to compute.They do not affect the asymptotic expansion of ϕ at r ≃ 1, ∞ but they increase the accuracy of the solution at finite r.Notice that all the arbitrary constants we are dropping upon integration correspond to trivial constant phase shifts.Plugging into equation 2.11 we get We are also in the position to check ϵ ≪ 1 given the explicit x we got.A direct computation reveals the condition which, as anticipated, puts a lower bound µ ≫ 2 3 √ 3 ≃ 0.4 on the l = 0 case.This reproduces the result obtained by Hui with a different technique.
How do we interpret this bound?In some sense, it comes from the significant smallness of the potential term µ 2 r 3 close to the horizon at small µ: the horizon is behaving similarly to a turning point 4 .In contrast, when l ≥ 1, the potential has developed zeros for µ even larger than 0.4, since µ c | l=1 ≃ 0.7: the breakdown takes place already at the µ c we computed in the previous subsection.We plot our approximate solution in figure 1.The accretion rate can be computed to be in agreement with ref. [44] The lower bound µ ≫ 0.4 compels us to study the complementary small µ regime µ ≪ 1, which we will tackle in the next subsection.
We anticipate that the connection problem will be solved by However, as made clear by figure 2, for l ≥ 2 there is a mass range not covered by either approximation.This will be dealt with in the appendix.
We anticipate that for l small, say l = 2, 3, the connection problem in this regime is approximately solved by

Small µr BH regime
Given the failure of WKB (and uniform approximations) in the small µr BH parameter range, we now look for a different kind of approximate solution to 2.3 for the case µr BH ≪ 1.The approach will be to neglect all terms proportional to µ in the differential equation and then identify the r domain where this is valid.We will then develop separate approximations where this assumption will fail.The upshot will be a subdivision of the domain in a near-horizon region, an intermediate region and a far region.
We refer to [44] for the l = 0 case, which agrees with our result.Our main contributions are the extension of the computation to all l and a more careful treatment of the matching of ϕ at different distances, which Hui et al. carried out simply by invoking continuity and derivability of the solution 'close' to the breakdown point.The technique we employ goes under the name of boundary layer theory (see the book by White [55]), but we will not assume any prior knowledge of the topic.The large mass approximation is valid in the blue area, while we can expand for µ ≪ 1 in the green area.Already at l ≥ 2 there exists a region beyond the scope of both approximations, which we cover in the appendix.This region enlarges if we don't want to approach the boundary of the blue and green areas, where our approximations become unreliable.
Since we work in the µ ≪ 1 approximation, we can expect µ 2 r 3 to be negligible in equation 2.3.This is true unless the r(r − 1) term coming from the angular momentum barrier becomes comparable or smaller.This happens in two regions, where we will instead approximate the r dependence: The attentive reader will recognize 1 + µ 2 and 1/µ 2 as the zeros of V (r) when µ ≪ 1.The issue we had in WKB with the zeros of V (r) has an avatar here as well.
Let us begin with the intermediate r region 1 + µ 2 < r < 1/µ 2 , where µ 2 r 3 is negligible.The solution can be written explicitly as where 2 F 1 is a hypergeometric function.
We now come to what we will call the two boundary layers: the near horizon region and the far region, defined by 2.23.In the near horizon region we can solve the approximate equation where we already imposed causality.Notice that we neglected the angular momentum part of the differential equation, but this is consistent if we don't keep track of (r − 1) −iµ+1 terms in the solution.Boundary layer theory also gives a rigorous criterion to systematically decide which terms one should keep, but it is just as easy to verify this a posteriori.
Finally, the far region is controlled by whose solution is where J and Y are Bessel functions.The reader should appreciate that the near horizon expansion keeps µ ln(r − 1) finite, while the far region is at finite µ √ r.In both cases there is a double scaling limit involving µ and r.Subleading terms are negligible with respect to these finite quantities.
Boundary layer theory also prescribes the (intuitive) matching condition: the near horizon layer is matched to the finite r solution under the assumptions µ ln(r − 1) ≫ 1, r − 1 ≪ 1.All terms that we can compute in both regions should coincide.
The near horizon and intermediate r expansions are where (a) l is the Pochhammer symbol.Correctly, all terms either match or are unknown (e.g. Similarly, we now match at large r (r ≫ 1, µ 2 r ≪ 1).We obtain the expansion of the scalar in the far region and in the intermediate region We now match the two by first obtaining some simple estimates: matching terms of order The conclusion is that we can safely forget about c 4 .We finally get Lastly, expanding the Bessel function at r → ∞, we obtain which has an overall µ scaling of µ −3/2−2l , thus growing with l.The physics underlying this result is simple: at higher angular momentum, light particles can rotate around the black hole, thus significantly enhancing the field.Indeed, this scaling can be deduced from a simpler argument: neglecting dimensionless factors, the leading-r behaviour in the intermediate-r regime is r l , which at the matching point r ∼ µ −2 amounts to µ −2l .Similarly, µ − 3 2 −2l /r 3/4 ∼ µ −2l at the same matching point.Due to analytic continuation in l, the result also agrees with the one obtained by Hui et al. for l = 0.

Love numbers
It is worth pointing out that, for intermediate distances, ϕ behaves as a massless scalar would, up to subleading µ corrections.The large r behaviour of ϕ intermed showcased in equation 2.33 suggests, for a massive scalar perturbation with µr BH ≪ 1, a definition of the Love numbers mirroring the one for the massless case: c 2 /c 1 (computed in eq.(2.31)).Since the result will turn out to be complex, this quantity will be called the scalar response coefficient (its real part being the actual Love number, the imaginary part being related to dissipation), and was calculated for example by Dubovsky et al [49] (their equation 6.17) for rotating black holes but a strictly massless scalar.Riotto et al. [50] were able to reproduce this result as well, using near horizon approximations (see their equation 5.24).We obtain which is also identical to their formula after some mathematical manipulations and after choosing ω = µ, as we do throughout.This agreement is not surprising, because the mass µ of the scalar only affects the solution at large distances (the frequency ω terms dominates the near horizon zone, while the mass µ dominates at large distances; the factor µ in our equations 2.31,2.36would more generally be ω).It makes sense that, in the limit µ → 0, the massless behaviour is recovered (at finite distances).The are however two nontrivial physical points: 1.When µr BH ≪ 1, there is an intermediate r region where the field behaves as if it were massless, allowing us to define the response coefficients and Love numbers.In the opposite regime µr BH ≳ 1, we do not have this parametric separation of scales and defining these quantities becomes much harder.We do not attempt to do so in this work.
2. Our method, besides its simplicity and generality, allows the computation of subleading corrections to 2.36 in ω and µ (even taken to be different).
We summarize the behaviour of ϕ in the three regions in table 1.We showcase our results by plotting the real and imaginary part of ϕ, this time for µ = 0.1.In order to have smooth transitions We plot in figure 3 the function which has all the correct asymptotic limits.An identical story plays out for the near horizon layer.
The accretion rate for the small mass regime was computed by Hui et al. in [44], who expressed it in terms of the DM density at intermediate distances for l = 0. Expressing it instead in terms of the density at large distances r ≫ µ −2 and dropping O(1) factors we get as a function of l where A is the amplitude of the wave at the horizon ϕ = A(r−1) −iµ .The rate is heavily suppressed as l increases, again due to the centrifugal barrier kicking in for small masses.

Numerical checks
We check the above formulas against the numerical results obtained in [51] by Hui et al.In their figure 9, the authors plot the inverse of the amplitude at r = 400r BH as a function of µ.They normalized the amplitude at the horizon to 1 and took 400r BH to be far enough that all our far  1) for µr BH = 0.1.We normalized the amplitude to 1 at the horizon.Due to the vast differences in the amplitudes, we plot each l separately.
region approximations should apply: indeed the largest length scale that appears in our problem is 1/µ 2 , which is much smaller than 400 for µ ≫ 0.05.
Their plot roughly captures the plateau at large µ, however in the small mass regime the ratio they plot goes like ∼ 1 cos(2µ √ r) | r=400 , meaning that one expects large fluctuations when numerically sampling different values of µ.To circumvent this problem, the authors actually average over small intervals in µ.This is a potential drawback of that plot, which ends up depending on the smoothing procedure.
Let us now move on to figure 10 of [51], where the authors plot |ϕ| 2 , normalized to 1 at the horizon.They fix l = m = 2, although m is irrelevant in the spherically symmetric case we are dealing with.The first plot is for µ = 0.01, meaning that the profile should accurately be described by our intermediate r region, which extends from r = 1 + 10 −4 to 10 4 .The analytic expression is Correctly, their plot presents no oscillations.The order of magnitude of the field is also correctly reproduced by this formula.
Next, let us look at the µ = 0.5 case.Here we are always in the large distance regime (r ≫ 24) and we expect our small mass approximation 0.5 ≪ 1 to give results correct up to O(1) corrections.We observe instead that the plot in [51] greatly underestimates, by a factor ∼ 10 4 , the field at large distances.To explain this discrepancy, observe that the authors seem to have normalized the field too far from the horizon, missing a sharp drop in the dark matter density close to r BH .
For µ = 1 the discussion in appendix A implies that intermediate masses behave roughly like large masses, hence our agreement with the decreasing profile.
Finally, let us look at figure 11 in [51], where they plot 1/|ϕ| at r = 400r BH for different values of µ.There are roughly three regimes: large mass, small mass and r ≫ 1 µ 2 , small mass and r ≪ 1 µ 2 .The authors provide estimates for all three regimes in their formula (4.1).Our formulas agree except when r is in the far regime r ≫ 1 µ 2 ≫ 1, which is their intermediate µ region, where the log-plot linearly grows.In our approximation |ϕ| should be given by so the slope 3 2 + 3l in [51] should be replaced by 3 2 + 2l, not as steep.A simple argument in favour of this is that at r ∼ µ −2 the last two formulas in 2.40 should parametrically coincide.This works well when µ is small, but the agreement with the numerical solution decreases as the graph plateaus on the right; there the approximation by Hui et al works better.We plot our approximation in figure 4.
The transition to the plateau on the right is predicted to happen when µ reaches µ c given by formula 2.7, thus the plateau region should move to the right as l increases.Contrary to what we see from the approximations in [51], the transition between the plateau on the left and the  [51].Referring to the equations for µ ≪ 1 in table 1, the green line is the intermediate r regime, the orange line is the far regime and the blue line is the very far regime once we set the cosine to 1, to tame the spikes.The transition between the green and the blue line is smooth if we consider the interpolating orange line.
region with linear growth in figure 4 is smooth once we properly take into account the large r regime (see table 1) described by a Bessel function.These examples demonstrate how employing a simple analytic approximation of the solution offers better control over the results compared to the challenges associated with using the full Heun function.

Conserved current
We pause to introduce a quantity we can prove is conserved on the equations of motion (very closely related to the wronskian).The beauty of the statement is that it is rigorous, exact, and it allows us to prove some claims.The downside is that it will not be enough to solve the problem.We define A simple substitution reveals dJ dr = 0 on-shell (2.42) Notice that J is conserved as a function of r (not time).Evaluated near the horizon on the causal (i.e.purely infalling) solution, one gets which, up to a factor, is the energy flux the wave carries into the black hole.
A simple thing we can check is that, if ω < µ, bounded solutions are not allowed.WKB (which is locally accurate away from turning points) gives ϕ ∼ e − √ µ 2 −ω 2 r at large distances.For this solution, J equals 0 due to the exponential suppression.Since the causal solution has J ̸ = 0, it never matches with a bounded profile for energies lower than the mass.Physically, this translates in the absence of (stable) bound states.We will later devote much discussion to solitons, which are bound states that evade this argument by having complex ω.
We can additionally check whether our approximate solutions respect this conservation law.The answer is trivially yes for the large mass regime approximation.However, in the small mass regime there seems to be a contradiction because J = 0 at large distances.Is our solution inaccurate?
The answer lies in the size of the ingoing and outgoing wave at large distances: within our accuracy the two amplitudes are large (∝ µ −3/2−2l ), which would naively suggest an equally large J at infinity.However because the amplitudes are precisely equal at this order of approximation, the waves sum to a cosine and the leading contribution to J cancels.We expect subleading orders in ϕ to lift this degeneracy and give a value of J much smaller than the individual contributions, consistent with the conservation of the current.We conclude that our approximate solution correctly captures the leading order behaviour of the exact solution, passing a nontrivial consistency check.
We will now see that J allows us to compute this slight difference in the amplitudes far away, without the need of computing the aforementioned subleading orders of ϕ.Very generally, where the equal sign is an exact statement, while the ≃ indicates that we are neglecting higher order terms which however do not affect J. 35.Then simple algebra gives

Imposing conservation one obtains
Therefore, as promised, J allows us to compute the very small difference |A − | − |A + |.
Dark matter solitons are not topologically protected: we will use the word to loosely mean any stationary field configuration in a bound state, thus the density of dark matter in a compact region will be significantly greater than in the rest of the halo.We will assume the following hierarchy of scales: The scenario we will consider is that the supermassive black hole is at the centre of the soliton.Given this hypothesis, the behaviour of such solitons can be dominated by the gravity of the supermassive black hole or by the self-gravity of dark matter.We begin the discussion from the first, simpler case.
What is the difference between the discussion of black hole domination of the previous section and the following, which we will call black hole dominated soliton?Self-gravity is still negligible but now we will consider bound states, implying "ω < µ ".Actually ω will turn out to have a very small imaginary part, evading the argument in section 2.4.We will first neglect the imaginary part and later compute it and check this assumption.
Denoting ω = µ(1 + γ) and expanding in small γ, the potential in equation 2.3 now has a term ∝ γr 4 .We will assume this term to be small enough that only the tail of the solutions found in the previous sections are modified, at a distance r ≳ |γ| −1 .Seeking only the behaviour of the solution at large distances, we can rewrite 2.3 approximately far from the horizon as To emphasize the scale r = |γ| −1 , we define r ≡ x −γ , obtaining This equation can be solved exactly, but a better intuition for the different distance regimes comes from a WKB approximation.To take this route, we now assume µ 2 |γ| ≫ 1.Given |γ| ≪ 1, this condition is obvious unless µ is in the small mass regime.The condition is equivalent to 1 |γ| ≫ 1 µ 2 , where 1 µ 2 is the length scale at which the decaying r −3/4 behaviour starts, in agreement with our hypothesis that only the tail of the solution should be affected.We shall further comment on this hypothesis later.Such a large parameter invites a WKB approximation, which will be accurate far from the turning point x = 1  2 , where the Airy matching will be employed.Then the decaying solution at x > 1/2 will be and at x ≪ 1/2 it matches We recognize the r −3/4 and the oscillating behaviour we already saw in the small mass regime (equation 2.35).The main observation is that if µr BH ≫ 1 there is a surprising clash between causal boundary conditions at the horizon, which demand oscillations that go like e −2iµ √ r (equation 2.14), and boundedness at infinity, which imposes a cosine (equation 2.51).Black hole dominated solitons can thus only be formed in the small µr BH regime.
We give the explicit very large distance behaviour of ϕ, normalizing the amplitude at the horizon to 1 as usual (that is, matching 2.51 to the very large r regime in table 1) We now estimate Im(γ) and check that it is small.The idea is to estimate black hole accretion, which is linked to the decay of the soliton.Neglecting O(1) factors so, using GM BH r sol µ 2 ∼ O(1), Finally we determine whether our hypothesis (µr BH ) 2 ≫ |γ| was well motivated.Finding γ amounts to determining the spectrum of bound states.To this end relying on an exact solution of equation 2.49 is more appropriate.Demanding boundedness at x = 0, vanishing boundary condition at x → ∞ and writing down the full solution in terms of 1 F 1 and U hypergeometric functions, we get the quantization condition suggesting that γ should indeed always be small.Had we not known about the exact solution, the quantization condition coming from WKB would have given us µ 2 |γ| = 8(n + 1 4 ) 2 .Given this good agreement, in appendix C we will trust the WKB result when computing the soliton spectrum in self-gravity domination, which is what we turn to in the next section.

Including self-gravity
A more interesting scenario is instead the soliton dominated regime, meaning that self-gravity dominates over the black hole gravity.We expect this regime to be roughly characterized by the condition M BH ≪ M sol , where we define the total soliton mass to be M sol = ∞ rBH ρ DM (r)4πr 2 dr (convergence is guaranteed by the exponential decay at infinity).Indeed we will see that black hole gravity is not necessary to hold together the soliton: self gravity of dark matter can be enough and we could meaningfully study the case M BH = 0.
One could think of the following possibility: a black hole surrounded by a soliton so massive that the dynamics is entirely dominated by dark matter self-gravity at all distances outside the horizon.We will show that this is actually impossible (assuming stability of the overall system), and that there is always a region outside the horizon dominated by the black hole and general relativity, however small M BH may be.

Scaling argument for the size of the soliton
We now include the self gravity of the scalar as a (non-relativistic) newtonian potential, which will obey a separate Poisson equation.The aim of the following discussion is to find an estimate for the soliton size.This discussion is fairly standard in the literature, and we will follow [14,32].Once we include self gravity, negative energy solutions appear (bound states).Energy will be minimized by a spherically symmetric field configuration, so we set l = 0 from here on.
Since it would be confusing to have both the scalar field ϕ and the newtonian potential Φ around, we define χ as The √ 8πG factor is chosen to simplify Poisson's equation.Notice how the energy density would be halved if the scalar were real.From now on we will always keep r BH explicit.Equation 2.3 is then amended to where by assumption Φ vanishes at large distances and ω = µ(1 + γ).We emphasize that this equation retains a fully relativistic treatment of the black hole in terms of its background metric, on top of which a correction is introduced to account for the self-gravity of the scalar.For a non-relativistic stress-energy tensor, this metric correction is suppressed by powers of c except for the time-time component, which is identified as the Newtonian potential.We further assume that Φ only affects the dark matter profile far from the black hole, so that the background metric in Poisson's equation can be approximated as flat.We will later observe the self-consistency of this assumption.At large r equation 3.2 reduces to which is the form used in [32].
In the approximation where the black hole is negligible, a well known scaling symmetry allows one to set χ(0) = 1.Then, numerical simulations [31,56,57] indicate that the smallest possible γ is ≃ −0.69 which partially supports the approximation |γ| ≪ 1. Scaled solutions turn out to have larger, less massive solitons and smaller γ.
To obtain the soliton size, we argue that GM sol ∼ Gρ s r 3 ∼ µ 2 χ 2 r 3 , while from 3.3,3.4we get Φ ∼ µ 2 r 2 χ 2 , χ ∼ µ 2 r 2 Φ χ .Substitutions yield GM sol r sol µ 2 ∼ 1, which can be compared with the analogous result obtained in black hole domination A much more detailed discussion of soliton domination can be found in appendix C. We notice that for M sol ≃ M BH the two formulas for r sol parametrically coincide.

Exact results
The purpose of this section is twofold.The first will be to exploit some mathematical results already available in the literature and apply them to a simplified version of our problem.Our aim will be to re-derive the same results using our WKB approximation, which will give us confidence in the results we will later find for the full problem under the same approximations.Secondly, we will derive an exact result (essentially the energy balance equation) which we will use to rigorously show that ω must be complex for bound states.This implies that bound states of the scalar will be slowly eaten away by the black hole, so they can only be meta-stable.
Therefore, we will now adapt some exact results and mathematical theorems, derived by Tod and Moroz in [52], to our problem.They will be far from enough to solve it, but they will enable us to put some of the our findings on firmer grounds (e.g. the quantization of the soliton spectrum).We discuss the simplified problem of a real scalar and no black hole.
As a first step, we map the system of differential equations 3.2,3.3onto a problem studied by Tod and Moroz in [52].They consider a system of differential equations for what they call S and V , which map to our variables through They assume S ∈ R (while we work with a complex field) and S(r) bounded and vanishing at infinity.They also assume r BH = 0, so there is no black hole.Their results can be summarized as follows: • There is a discrete family of solutions labelled by n ∈ N, the n − th solution having n − 1 zeros and vanishing at large r as • γ n < 0 and increasing with n towards zero • There are no other solutions To develop a better understanding of Tod and Moroz's findings, we begin by rederiving their results using our WKB approximation.
The WKB solution of 3.4 can be immediately written.Defining C 2 = 1 − (1 + γ) 2 , we have for A, B constants.We still don't know Φ, but supposing −2Φ − C 2 > 0 at large distances and recalling that Φ vanishes at infinity, we would get χ ∼ 1/r and the total mass would diverge.Thus −C 2 < 0 and γ < 0.
Thus −2Φ − C 2 < 0 at large r, and imposing boundedness we get which has the predicted asymptotic behaviour.As we decrease r, we can expect 2Φ + C 2 to change sign at some finite distance that we interpret as the soliton size r sol .
The solution at r < r sol will transition from exponentially decaying to oscillating and the coefficients A, B will (in principle) be calculable given A ′ and γ thanks to the WKB connection formulas.Generically 3.7 will diverge at the origin, unless A = −B.This condition quantizes the frequency spectrum, giving rise to the set of γ n .As we increase γ, C 2 decreases approaching zero, while r sol increases.In analogy with what happens in quantum mechanics, the number of zeros of the solution increases by one for each new bound state.This concludes the derivation of the results in [52] using WKB.Notice how we chose A ′ to be real.We can however rotate it by an overall phase, thus generating all bounded solutions of the complex scalar case with a simple overall phase rotation of χ.Having overcome this limitation, we now turn to the second one, namely the lack of a black hole which we now introduce back.
We now reconsider equations 3.2 and 3.3.They were obtained by treating the black hole within general relativity, and adding the scalar self gravity as a weak newtonian potential.We will now write down the equation for the scalar field and generalize the argument given in section 2.4 to prove that bound states are metastable: their frequency ω cannot be purely real.We do this by generalizing the previously conserved current: the main result is that J is not conserved anymore, and the deviation from dJ dr = 0 is proportional to the decay rate of the soliton.
Assuming spherical symmetry we can choose g θθ = r 2 , g ϕϕ = r 2 sin 2 θ as usual.We expect the scalar field ϕ to deform the Schwarzschild metric and lead generically to some g tt = −A(r), g rr = B(r).The equation for the scalar field is and J can be straightforwardly defined and shown not to be conserved We now want to show that ω cannot be real for soliton configurations.From equation 3.10 real ω implies conservation of J; however a direct computation reveals that bound states have J = 0 at infinity due to the exponential suppression of ϕ, while ϕ close to the horizon is still a purely infalling wave, hence J ̸ = 0. Since this violates dJ dr = 0, we have proved our claim Im(ω) ̸ = 0.
This fact is sometimes overlooked in the literature, however it has a clear physical counterpart: bound states of the scalar will be slowly eaten away by the black hole, so they cannot be stable.Actually, Im(ω) is the inverse time in which the bound state decays.
We leave it for appendix B to present and solve a concern one can have over the way we impose causality at the horizon.When the frequency is complex, ϕ causal diverges at the horizon while ϕ non-causal goes to zero.Since the non-causal solution is already exponentially suppressed at the horizon (and strictly zero at r = r BH ), one might worry that the causality condition has to be rediscussed.Actually, a more careful analysis using wave packets shows that the physical solution is always bounded and the usual causality condition should be imposed.This is reminiscent of black hole quasinormal modes, whose solutions are also divergent in frequency domain.

Solution from large to small r
We now explore the solution of 3.4 as we vary r.If the black hole dominates at all distances, we can neglect Φ when determining χ.As a second step, one can then determine Φ by solving 3.3.Since Φ does not affect the scalar profile in this scenario the problem is reduced to the one solved in section 2.
In this section we instead assume that the soliton dominates the dynamics at large distances, meaning M BH ≪ M sol .One goal will be to estimate Φ and consequently understand if the assumption r BH ≪ |rΦ| breaks down.It turns out that it always breaks down, at a distance we will name r e .At shorter distances we will then match the solution to the one where we assume r BH ≫ |rΦ|, which is the regime discussed in section 2. Finally, we assume Im(ω) to be negligible compared to µ, since its inverse is the time scale of the soliton decay, which we assume to be long-lived.We will check the self-consistency of this claim once we obtain the decay time of the soliton.
We start with the WKB solution of 3.4 for χ.Since the soliton dominates, M sol ≫ M BH and the results obtained by Tod and Moroz should approximately apply.The WKB approximation was able to reproduce all their results (see subsection 3.2), so we will trust this solution until the condition r BH ≪ |rΦ| is no longer met.
Recall that we labeled the size of the soliton r sol , which in WKB language coincides with the (single) zero of 2Φ + C 2 .The hierarchy of scales (which will be proved in the following) can be summarized as We can write the WKB ansatz for the solution both at r > r sol where it decays and at r < r sol where it oscillates.We then match them using the well known WKB connection formulas, resulting in , r < r sol (3.12)where C 2 = 1 − (1 + γ) 2 .We take this to be accurate away from the turning point, where the denominator (−2Φ − C 2 ) becomes zero and WKB breaks down while the true solution remains of the same order of magnitude.
Given our experience from the case without a black hole, if we consider the n-th bound state, the cosine oscillates n times.Besides, to restrict to physically meaningful solutions we should impose cos = 0 at the origin, so that χ remains finite.To restrict to the ground state, we should therefore be able to write with the understanding that µ We now make an important simplifying assumption (we will check its validity later): we assume −2Φ − C 2 to always be of the same order k when r e ≪ r ≪ r sol , so that χ ∝ sin(µkr) r and we can compute ρ(r) as , for r ≳ r sol (3.14)where O(1) factors (not essential) are determined from continuity and imposing that the total mass should indeed be M sol .The region r ≲ r e will be discussed later.Our approximation neglects the halo which is hosting the soliton.As a consequence, we are unable to discuss the relation between the halo and the soliton mass found in [12,29,31].However we expect the halo to have a negligible effect on the solution, as long as the soliton dominates over the halo density.
We are now in the position to carry out the program outlined at the beginning of this subsection: computing Φ and estimating r e .From Poisson's equation, Φ can be computed starting from infinity and imposing the boundary condition Φ(r) ∼ − GM sol r .It turns out that, as expected, the potential is well approximated by this expression everywhere outside the soliton, despite the exponential tail of the density.Inside the soliton We now determine c 1 , c 2 .Not only does c 1 = 0 make Φ ′ (r sol ) continuous (up to 8% error), but it makes Φ finite at the origin as we expect.c 2 ≃ −1 can be deduced from the continuity of Φ(r sol ).
A more drastic approximation would have been to average χ 2 ∼ cos 2 r 2 ∼ 1/2 r 2 , thus giving a potential which goes like Φ ≃ GM sol r sol (ln(r/r sol ) − 1) inside the soliton.This is in good agreement with the more complicated potential we got in equation 3.15, except close to the origin, where only the more accurate approximation shows that 2Φ + C 2 changes slowly and replacing it with a constant should yield reliable results.
We now determine the length scale r e , where −2rΦ ≃ r BH .Assuming M BH /M sol ≪ 1 we expect this scale to be much smaller than r sol .Exploiting this, we can expand rΦ close to the origin where it is very closely linear.We obtain Observe first that r e ≪ r sol and secondly that the black hole dominated region r BH < r < r e always exists, because stability of the soliton forces GM sol ≪ r sol which implies and therefore r BH ≪ r e .This confirms the whole hierarchy of scales we assumed in 3.11.
Is the stability requirement GM sol ≪ r sol trivially satisfied?As discussed in the review [14], it is not.In section 3.1 we linked the soliton size to the masses µ, M sol .This constrains the possible values of µ as follows: which therefore implies This is a very general bound that links the microscopic quantity µ with the masses of astrophysical objects.If a soliton is observed in a galaxy that hosts a SMBH with mass M BH , it implies an upper bound on the mass of the scalar: For completeness, let us point out that solitons dominated by a black hole also have a stability bound µr BH ≪ MBH M sol , which however is trivially satisfied because, due to the discussion in 2.5, µr BH ≪ 1 in this regime.A novelty that we bring is that the bound for soliton domination also forces us to restrict to the small µr BH regime, where we know that boundary conditions imposed at the horizon become less and less important the further away we are from r BH (in fact the part of the solution that depends on them decays like ∼ 1/r compared to the part which does not).Given this and r e ≫ r BH , we deduce that boundary conditions are unimportant in the regime of soliton domination.Quantitatively, suppression is ∝ rBH re = GM sol r sol .

Solitons at small r
We now match the solution obtained thus far with what we got in section 2 under the approximation that the black hole dominates the dynamics, which happens when r < r e .From the bound µr BH ≪ MBH M sol ≪ 1 we can focus on the small mass regime only.
Far from the horizon the field can be approximated as We now show that µ √ r BH r e ≪ 1, meaning that we are close to maximum of χ.In the language of section 2, this is the intermediate r region.Some simple algebra and the definitions of r e and r sol lead to We can then approximate χ with a constant.The density is matched at r = r e , ρ(r) = ρ(r e ) = (8.3) 3.21 was derived assuming γ = 0, but this doesn't matter.In fact, going through the derivation of the wave profile one has to replace µ → ω in the expression of the wave near the horizon, which gives a negligible contribution if r ≪ 1 µ 2 rBH .Up to very small corrections ∝ , this density is what one would obtain had they neglected the black hole altogether, which is the usual assumption in the literature.The point of our analysis was to check that the part of the solution dependent on the boundary conditions never enters a nonlinear regime.
We can thus fill in the missing piece of equation 3.14: We present the plot of this function in figure 5.
We end this section computing the accretion rate of the black hole.This quantity depends on the precise value of ρ at the horizon.Matching ρ hor to the soliton profile outside the sphere of influence of the black hole was always achieved under some simplifying assumptions (i.e.absence of features in the dark matter density profile inside the black hole dominated region) we now have under control.Employing 3.24 we have We can rewrite r sol using GM sol r sol µ 2 ≃ O(1), where the numerical constant is computed in appendix C to be ≃ 6 (note that some authors report roughly half).Then we have which agrees with the numerical results in equation (11) in ref. [42] (their numerical prefactor corresponds to 0.06).
As promised at the beginning of section 3.3, we now check that Im(ω) ≪ µ.From 3.25 we can extract To see that this is ≪ 1, we square and note the following inequalities from the stability of the soliton.

Discussion
We reviewed black holes endowed with nonrotating wave dark matter hair and extended the analysis to l > 0. When µr BH ≪ 1 interference patterns typical of wave dark matter emerge, giving rise to O(1) fluctuations in the density profile.However, in the opposite regime, these effects are washed out by the black hole gravity.This last conclusion crucially relies on imposing causal boundary conditions at the horizon, suggesting that caution should be used when modelling a black hole using newtonian gravity, even at large distances.In addition to providing control over all distances, our techniques offer a precise method for calculating the dynamical Love numbers for scalar perturbations to leading order in small ω and µ, which can be easily extended to subleading orders.
Our primary focus is on dark matter solitons centered around a supermassive black hole.In a black hole dominated soliton, M BH ≫ M sol , the scalar profile is given by the solution presented by Hui et al. up to distances r ∼ r BH /|γ| beyond which it exponentially decays.As we discuss in section 2.5, such solution exists only in the small mass regime µr BH ≪ 1. Solutions with larger mass are incompatible with causal boundary condition at the black hole horizon.
In the opposite regime, M BH ≪ M sol , we analytically computed the approximate soliton profile at large (small) distances where gravity is dominated by the soliton (black hole), respectively.At the radius r e , representing the sphere of influence of the black hole, we matched the density profiles.The soliton's stability in this scenario enforces a hierarchy of scales, r BH ≪ r e ≪ r sol , which also implies a small mass limit, µr BH ≪ M BH /M sol .This condition indicates that causal boundary conditions at the horizon have no effect far from BH, similar to the small mass regime in black hole domination.
We emphasize that having analytical control over the solution enables us to perform these assessments.With a solution that remains valid even at horizon scales, we can reliably compute the accretion rate of the black hole.
We now address some limitations of our work.Firstly, we do not model the dark matter halo, which is a notoriously challenging task.However, we expect that our results remain applicable as long as the dark matter density within the soliton overwhelms the halo density.Performing a matching of our solution with the halo profile is beyond the scope of this work.
Secondly, we have focused on non-rotating black holes.We do not believe this assumption to be overly restrictive, given that Sagittarius is a slowly spinning black hole with J/M < 10%.In addition [51] numerically studied the issue and found moderate corrections for J/M up to 50%.Nevertheless, accounting for a (slowly) spinning black hole would provide an interesting cross-check of those results.Because the analogue of equation 2.3 for spinning black holes can be similarly written, we expect the analysis to be feasible using techniques akin to the ones applied in this paper.However, we leave this investigation for future research.
Lastly we neglected any non-gravitational interaction.A discussion of the validity of this assumption can be found for example in [32], where the authors conclude that a wide range of ALPs can satisfy this approximation.

Aknowledgments
We thank Lam Hui, Luca Santoni, Guanhao Sun, Giovanni Maria Tomaselli, and Rodrigo Vicente for helpful discussions and comments on the draft.

A Intermediate mass µ at l ≥ 1
This appendix is dedicated to developing an approximate solution of 2.3 for the white are in figure 2.
The small µ approximation is by definition invalid, so we would like to again resort to WKB.The main complication is that now V (r) has two zeros.To tackle this problem, we fully exploit the technique of uniform approximations we introduced in section 2.
We will first find a differential equation (that we can explicitly solve) with the same structure of zeros as 2.3 in the region r > 1.We will then map this exact solution to an approximate solution of the original 2.3.
We keep the definitions of r, x given in 2.4,2.10.We now also have the zeros of V (r), given by r 1,2 as defined in 2.6.
As auxiliary differential equation, we choose 2.9 where Γ should now capture the zeros of V (r).The simplest choice is Γ(x) ≡ k(x − r 1 )(x − r 2 ), k > 0. The zeros will truly match only if we impose (see [48] for a detailed discussion of this point) which we use to fix k.We postpone the issue of writing down an explicit formula for k.
We now solve 2.9 with the new choice of Γ(x) and we expand the solution in the near horizon and large distance regimes.Equations 2.11, 2.12 and 2.13 go through.We replicate the estimates for x at large r and close to the horizon, obtaining

B Causality at the horizon
This appendix is dedicated to a more careful analysis of the causality condition when ω is not real.We will assume stability, so that the imaginary part of ω must be negative Im(ω) < 0 since ϕ ∼ e −iωt and we want a decaying solution.
Now as usual with complex frequencies, ϕ diverges near the horizon.Indeed, the approximate solution is ϕ ≃ (r − 1) −iω ∼ (r − 1) −iℜ(ω) (r − 1) Im(ω) (B.1) while the non-causal solution is zero at the horizon ϕ ≃ (r − 1) iω ∼ (r − 1) iℜ(ω) (r − 1) −Im(ω) (B.2) Hui and Baumann impose that the amplitude of the non-causal solution should be zero at the horizon.Since in this setup this is so automatically, does this mean that we can forget the causality condition?On the same note, should we be worried that the amplitude of the causal solution diverges at the horizon?
Let us answer the question by schematically building a wave-packet out of the non-causal wave components B.2 (the story for B.1 is similar).We track the evolution of this (outgoing) wave packet in time and we ask what happens at very early times.If the wave packet vanishes close to the horizon causality is safe, otherwise it isn't.
We anticipate that the wave packet amplitude (at its peak) is almost unchanged for all t.This is so despite its wave components B.2 approaching zero for r → r BH at any given time t.Thus we want to set the non-causal solution to zero, otherwise an observer hovering outside the horizon would see stuff flying out.A free falling observer would instead perceive this as an infinite energy solution.Now we carry out the analysis sketched above.The wavepacket is a superposition of B.2 with modulation φ(ω).We are neglecting the infrared (large r) deformations of B.2, but they are irrelevant for the present purpose.ϕ(r, t) = dω e iω ln(r−rBH) e −iωt φ(ω) (B.3) As an example, we choose φ(ω) = e −iω ln r0 f ω0 (ω).The exponential prefactor translates the solution at time t = 0. f is some real function peaked at ω 0 , which gives the size of the wave packet.
We split real and imaginary part of ω into ω → ω − iϵ ϕ(r, t) = e ϵ ln(r−rBH)−ϵt−ϵ ln r0 dω e iω ln(r−rBH) e −iωt e −iω ln r0 f ω0 (ω) (B.4) and we ask when is the integral large.A saddle point approximation yields the position of the peak of the wavefunction in time: The wavepacket flies out of the horizon as predicted, and its initial position is controlled by r 0 as initially claimed.Let's track the wave packet (peak) amplitude in time.It is defined as A(t) ≡ ϕ(r(t), t) ∼ e ϵ ln(r(t)−rBH)−ϵt−ϵ ln r0 = 1 (B.7) so it remains finite (even constant) at all times.While constancy is an artefact of the near horizon limit (i.e. the amplitude does change at large t), finiteness is robust and implies that the wavepacket is non-zero arbitrarily close to the horizon.
Similarly, a wave packet build with the causal waves B.1 will also remain finite as it approaches the horizon, despite the individual waves blowing up at r = r BH . in order to have χ bounded at the origin and matched to an asymptotically decaying solution at infinity.This means that The first conclusion is that the approximate relation µ 2 GM sol r sol ≃ 1 found in 3.5 using a scaling argument is essentially correct, but it can now be made more precise.Reinserting ℏ, c we obtain for the lowest energy level n = 0 GM sol r sol µ 2 ≃ (6.2), M sol r sol ≃ (5.3) which is roughly twice the value reported in [32,34].The discrepancy could be due to the approximations we have made.Moreover, we have a prediction for the mass scale of all excited bound states n > 0 Additionally, Bar [32] has a prediction for what γ, M sol , r sol should be if we set χ = 1 at the origin (so that the scaling symmetry is broken).We will now try to reproduce their numerical results.
First, M sol and r sol are related by the already mentioned relation GM sol r sol µ 2 ≃ 1, so we will only determine r sol .We start from χ(0) ≡ χ 0 = 1, which implies Bar has µr sol ≃ 1.3 in [32], so again we are qualitatively correct.It also makes sense that for ground state of the soliton to have size roughly a few compton wavelengths of the scalar particle.The scale parameter λ, which we can define as χ 0 = λ 2 , thus has the further interpretation that λ −1 counts the number of oscillations of the field inside the soliton.
Finally, we now turn to the computation of C 2 .From C 2 = 2GM sol r sol , we can deduce where we first used our own µr sol ≃ 2.7 and then we used Bar's µr sol ≃ 1.3.To compare with the value reported in [31,56,57], we should make the (improper) identification C 2 ≡ −2γ, which brings our result to −γ ≃ 1.7/2 = 0.85, not far from 0.69 obtained numerically.

Figure 1 :
Figure 1: Plot of our large mass approximation of ϕ 2.18 for different values of µr BH and l.The amplitude is normalized to 1 at the horizon.

Figure 2 :
Figure 2:The large mass approximation is valid in the blue area, while we can expand for µ ≪ 1 in the green area.Already at l ≥ 2 there exists a region beyond the scope of both approximations, which we cover in the appendix.This region enlarges if we don't want to approach the boundary of the blue and green areas, where our approximations become unreliable.

Figure 3 :
Figure 3: Plot of ϕ (table1) for µr BH = 0.1.We normalized the amplitude to 1 at the horizon.Due to the vast differences in the amplitudes, we plot each l separately.

Figure 4 :
Figure 4: Plot of |ϕ(r = 1)|/|ϕ(r = 400)| as µ changes.l = 1 on the left, l = 2 on the right.This plot corresponds to figure 11 in[51].Referring to the equations for µ ≪ 1 in table 1, the green line is the intermediate r regime, the orange line is the far regime and the blue line is the very far regime once we set the cosine to 1, to tame the spikes.The transition between the green and the blue line is smooth if we consider the interpolating orange line.

,
for r e ≲ r ≲ r sol− GM sol r ,for r ≳ r sol(3.15)

4 )
Substituting the value of Φ obtained in 3.15 we can do the integral and obtain µ GM sol r sol (0

Table 1 :
Behaviour of ϕ in the small mass regime.