# Quasistationary solutions of scalar fields around accreting black holes

###### Abstract

Massive scalar fields can form long-lived configurations around black holes. These configurations, dubbed quasi-bound states, have been studied both in the linear and nonlinear regimes. In this paper we show that quasi-bound states can form in a dynamical scenario in which the mass of the black hole grows significantly due to the capture of infalling matter. We solve the Klein-Gordon equation numerically in spherical symmetry, mimicking the evolution of the spacetime through a sequence of analytic Schwarzschild black hole solutions of increasing mass. It is found that the frequency of oscillation of the quasi-bound states decreases as the mass of the black hole increases. In addition, accretion leads to a significative increase of the exponential decay of the scalar field energy due to the presence of terms of order higher than linear in the exponent. We compare the black hole mass growth rates used in our study with estimates from observational surveys and extrapolate our results to values of the scalar field masses consistent with models that propose scalar fields as dark matter in the universe. We show that even for unrealistically large mass accretion rates, quasi-bound states around accreting black holes can survive for cosmological timescales. Our results provide further support to the intriguing possibility of the existence of dark matter halos based on (ultra-light) scalar fields surrounding supermassive black holes in galactic centers.

###### pacs:

95.30.Sf, 04.70.Bw, 04.25.dgJuly 15, 2019

## I Introduction

There is compelling evidence that most nearby galaxies host supermassive black holes (SMBHs) in their centers, with masses in the range Ferrarese and Ford (2005); Kormendy and Ho (2013). Such SMBHs are expected to be the evolutionary result of the growth of seed BHs in high redshift galaxies through accretion episodes and mergers of massive BHs (MBH) binaries with masses somewhere in between those of stellar-origin BHs and SMBHs Kauffmann and Haehnelt (2000); Volonteri et al. (2003). The discovery of supermassive luminous quasars at redshifts of has shown that SMBHs with masses , must form extremely early on in the history of the universe Barth et al. (2003); Willott et al. (2005); Mortlock et al. (2011). These SMBHs must grow rapidly in order to acquire its mass within a short period of Gyr. Explaining the formation and evolution of SMBHs dwelling in the most powerful quasars when the Universe was less than 1 Gyr old (and of the regular and much smaller MBHs hidden in 13 Gyr old galaxies) is a pressing open issue Volonteri (2010). Proposed models involve the gravitational collapse of gas clouds or the collapse of supermassive stars in the early universe (a model hampered by the low masses of the initial seeds of first-generation Pop III stars Madau and Rees (2001)), the runaway growth by accretion onto Pop III BH, or mergers of smaller size BHs. In either scenario, the formation of SMBHs is a highly dynamical event amenable to gravitational wave investigations. Indeed, it is expected that eLISA will probe MBH binaries in the range out to redshift through the detection of their gravitational waves in the mHz frequency band Barausse et al. (2015).

While one can make a convincing case that SMBHs have grown largely through accretion, with the consequent energy emission observed in electromagnetic output, it has been argued that an exponential growth at the Eddington-limited e-folding time is too slow to grow stellar-mass BH seeds into the supermassive luminous quasars that are observed at Mortlock et al. (2011). Some proposals to circumvent this issue invoke super-Eddington accretion for brief periods of time Alexander and Natarajan (2014), the formation through accretion and merging of the first stellar remnants Bromm and Loeb (2003); Haiman and Loeb (2001) or via more massive seeds from the collapse of pre-galactic disks at high redshifts Begelman et al. (2006); Begelman (2010); Mayer et al. (2010).

In light of these proposals it is worth considering if rapid BH accretion may have any effect on the distribution of the associated dark matter halo when the latter is modeled as a scalar field. In the absence of accretion the existence of long-lasting scalar field configurations surrounding a (non-rotating) BH has been investigated in a number of recent papers, either in the test-field limit Barranco et al. (2011, 2012, 2014) or employing self-gravitating scalar fields Sanchis-Gual et al. (2015a, b). These papers have shown that SMBHs at galactic centers do not represent a serious threat to dark matter models based on (ultra-light) scalar fields as a viable alternative to the usual description of dark matter in terms of weakly interacting massive particles. For both, scalar fields around SMBHs and axion-like scalar fields around primordial BHs, it has been found that scalar fields can survive for cosmological timescales Barranco et al. (2012). Despite these findings, and due to the rapid growth of accretion-powered SMBHs, it is worth investigating the chances of survival of the scalar field within such dynamical situation. This is the purpose of this paper. Here we study the properties of scalar field quasi-bound states in the background of an accreting spherically symmetric BH. Assuming that the BH mass grows adiabatically due to infalling matter we show that indeed, long-lasting, quasi-bound states can still form in such scenario.

This paper is organized as follows: in Section II we lay out the mathematical and physical approach we use to carry out our investigation. In particular Section II.5 contains a brief synopsis of our numerical methodology. The results are presented and discussed in Section III, while Section IV contains the summary of our findings. Throughout the paper we employ geometrized units (). Latin indices indicate spatial indices and hence run from 1 to 3 while Greek indices run from 0 to 3.

## Ii Problem setup

### ii.1 Klein-Gordon equation

Our setup considers a scalar field distribution of sufficiently small energy to neglect its self-gravity, i.e. the field can be regarded as a test-field. This configuration surrounds a BH which is assumed to be continuously accreting matter. The dynamics of the scalar field is governed by the Klein-Gordon equation,

(1) |

where the D’Alambertian operator is defined by . We follow the convention that is dimensionless and , the mass of the scalar field, has dimensions of (length).

We write the spacetime metric as

(2) | |||||

where is the lapse function, the shift vector, and the spatial metric. We adopt a conformal decomposition of the spatial metric

(3) |

where is the conformal factor, the conformally related metric and its determinant. By assuming spherical symmetry the line element may be written as

(4) |

with being the solid angle element and and two independent metric functions.

To solve the Klein-Gordon equation in spherical symmetry, we introduce two first-order fields, and , defined as:

(5) | |||||

(6) |

where is the unit vector normal to the surfaces of constant . Therefore, using Eq. (1) we obtain the following system of first-order equations:

(7) | |||||

(8) | |||||

(9) | |||||

where is the trace of the extrinsic curvature. The stress-energy tensor of the scalar field reads

(10) |

From this tensor we can compute the energy of the scalar field

(11) |

where is the radius of the apparent horizon and is given by

(12) |

### ii.2 Analytic black hole solution

The numerical evolution of a single BH using the so-called “moving puncture” technique leads to a well-known time-independent, maximal slicing solution of Schwarzschild Hannam et al. (2007). It was shown in Baumgarte and Naculich (2007) that this solution can also be constructed analytically and be used as a test for numerical codes. We take advantage of this result to put forward the defining characteristic of the procedure we employ in the current investigation. Namely, we avoid evolving numerically the BH, computing instead a sequence of analytical solutions at each time step for different BH masses to mimic the BH growth for different accretion rates. This procedure allows us to simulate long episodes of accretion without resorting to a test fluid in order to achieve high accretion rates.

The analytic solution is constructed by defining all quantities as a function of the gauge-invariant areal radius . We have to transform the solution into isotropic coordinates, as the latter are used in our numerical procedure, comparing the spatial metrics as a function of and the isotropic radius

(13) |

Thus, the isotropic radius is given as a function of by

(14) | |||||

In the limiting case , . Correspondingly, the conformal factor is obtained from

(15) | |||||

Finally, the lapse function and the isotropic shift vector are respectively given by

(16) | |||||

(17) |

### ii.3 Adiabatic growth of the BH mass

We assume that the mass of the BH grows due to the capture of matter as e.g. falling in from an accretion disk. This infalling matter is assumed to interact with the scalar field only gravitationally, that is the quasi-bound states surrounding the BH can only be affected by the increase of the mass of the BH. We will employ a simple phenomenological law based on observational and theoretical grounds which allows us to reasonably incorporate the growth of the BH mass in our model. As mentioned before, the observational evidence of the existence of very luminous quasars, which implies BH masses of at Willott et al. (2010); Fan et al. (2006); Mortlock et al. (2011) demonstrates that SMBHs grow rapidly in a short span of time ( years). Moreover, cosmological simulations Springel et al. (2005); Hu et al. (2006); Di Matteo et al. (2008) suggest that SMBH seeds undergo an exponential growth phase at early times, . Therefore, given a growth rate , the mass of the BH will increase as

(18) |

where is the initial BH mass and is the time in our simulations. The actual mechanism that produces the growth of the BH mass is not actually relevant for our study since we are interested in describing the evolution of the scalar field. Thus, assuming that it grows according to Eq. (18) seems quite convenient.

We first consider five different values for , namely, and . The time evolution of the BH mass for the different growth (accretion) rates considered is plotted in Fig. 1. We note that our values of , while small, are nevertheless orders of magnitude higher than realistic values inferred from observations, which are at most in our units Tsai et al. (2015). (We further expound on this issue in Section III.3 below.) Using such values would however render the numerical investigation prohibitively expensive. The first set of simulations reported in this paper evolve the scalar fields up to a final time , which is s for . Despite the evolutions are fairly long from the computational point of view, they are certainly not so on astrophysical grounds. Therefore, in order to study the effect of accretion on the scalar field evolution in affordable computational times, we have to resort to large enough growth rates. As we show below, our simulations indicate that even for the smaller considered some slight differences appear by the end of the simulation. Nevertheless, longer, and computationally more expensive, simulations would be necessary to show the influence of the growth rate for the smaller values of .

In the second part of our study we change the computational setup of the problem by placing reflecting boundary conditions for the scalar field at some radius, i.e. placing a mirror beyond which the scalar field is required to vanish, as we did in Sanchis-Gual et al. (2016). While evolving the scalar field in a cavity is certainly an unrealistic situation, it has some practical advantages as it allows us to investigate lower values of and perform longer runs since the computational grid is significantly smaller than in the “natural” setup (with outgoing boundary conditions at the outermost radial grid zone). Within this idealized setup the initial value of can be as low as . Despite this reduction only brings in practice about one order of magnitude gain in the final evolution time, it is nevertheless significant to better quantify the exponential decay in the energy of the scalar field, as we show below.

### ii.4 Initial data

As initial data for the scalar field we set a Gaussian distribution of the form

(19) |

where is the initial amplitude, is the center of the Gaussian, and is its width. The auxiliary first order quantities are initialized as follows

(20) | |||||

(21) |

Likewise, we choose a conformally flat metric with together with a time symmetry condition, i.e. a vanishing extrinsic curvature, .

For the BH, we compute the Schwarzschild solution using the moving puncture technique from equations (14)-(17) for the initial mass .

1 | 2 | |||||||
---|---|---|---|---|---|---|---|---|

0.00 | 0.05 | 9.42E-05 | 0.04998 | … | 3.257E-08 | 0 | 0 | 1.00 |

5E-07 | 0.05 | 9.42E-05 | 0.04993 | … | 3.272E-08 | … | … | 1.02 |

5E-06 | 0.05 | 9.42E-05 | 0.04989 | … | 2.390E-08 | 2.263E-13 | … | 1.22 |

5E-05 | 0.05 | 9.42E-05 | 0.04914 | 0.04951 | 4.426E-06 | -4.755E-10 | 1.548E-14 | 7.39 |

1E-04 | 0.05 | 9.42E-05 | 0.04844 | 0.04908 | 9.826E-06 | -1.787E-09 | 9.803E-14 | 54.60 |

0.00 | 0.08 | 6.18E-05 | 0.07969 | … | 2.179E-06 | 0 | 0 | 1.00 |

5E-07 | 0.08 | 6.18E-05 | 0.07969 | … | 2.284E-06 | … | … | 1.02 |

5E-06 | 0.08 | 6.18E-05 | 0.07959 | … | 1.886E-06 | 4.773E-11 | … | 1.22 |

5E-05 | 0.08 | 6.18E-05 | 0.07905 | 0.07942 | 2.571E-05 | -4.547E-09 | 3.318E-13 | 7.39 |

1E-04 | 0.08 | 6.18E-05 | 0.07825 | 0.07895 | 5.309E-05 | -1.541E-08 | 1.693E-12 | 54.60 |

0.00 | 0.10 | 5.00E-05 | 0.09945 | 0.09989 | 1.494E-05 | 0 | 0 | 1.00 |

5E-07 | 0.10 | 5.00E-05 | 0.09944 | 0.09988 | 1.551E-05 | … | … | 1.02 |

5E-06 | 0.10 | 5.00E-05 | 0.09928 | 0.09978 | 1.935E-05 | 7.147E-11 | … | 1.22 |

5E-05 | 0.10 | 5.00E-05 | 0.09878 | 0.09929 | 6.962E-05 | -1.255E-08 | 1.455E-12 | 7.39 |

1E-04 | 0.10 | 5.00E-05 | 0.09808 | 0.09873 | 1.054E-04 | -2.974E-08 | 5.194E-12 | 54.60 |

0.00 | 0.15 | 3.37E-05 | 0.14801 | 0.14950 | 3.394E-04 | 0 | 0 | 1.00 |

5E-07 | 0.15 | 3.37E-05 | 0.14793 | 0.14954 | 3.399E-04 | … | … | 1.02 |

5E-06 | 0.15 | 3.37E-05 | 0.14778 | 0.14943 | 3.450E-04 | … | … | 1.22 |

5E-05 | 0.15 | 3.37E-05 | 0.14785 | 0.14892 | 3.950E-04 | … | … | 7.39 |

1E-04 | 0.15 | 3.37E-05 | 0.14745 | 0.14807 | 4.360E-04 | … | … | 54.60 |

0.00 | 0.20 | 2.54E-5 | 0.19879 | 0.19948 | 6.800E-04 | 0 | 0 | 1.00 |

5E-07 | 0.20 | 2.54E-5 | 0.19882 | 0.19955 | 6.800E-04 | … | … | 1.02 |

5E-06 | 0.20 | 2.54E-5 | 0.19868 | 0.19943 | 6.830E-04 | … | … | 1.22 |

5E-05 | 0.20 | 2.54E-5 | 0.19830 | 0.19855 | 7.000E-04 | … | … | 7.39 |

1E-04 | 0.20 | 2.54E-5 | 0.19707 | 0.19787 | 7.197E-04 | … | … | 54.60 |

### ii.5 Numerical approach

The solution of the Klein-Gordon evolution equation is computed with the same type of numerical techniques we have used in previous works. The reader is addressed in particular to Refs. Montero and Cordero-Carrion (2012); Sanchis-Gual et al. (2015a) for full details on those techniques. As a succinct summary we mention that the evolution equations are integrated using the second-order PIRK method introduced in Cordero-Carrión and Cerdá-Durán (2012, 2014) which allows to handle singular terms that appear due to our choice of curvilinear coordinates. The spatial derivatives in the evolution equations are computed using a fourth-order centered finite difference approximation on a logarithmic grid except for the advection terms for which we adopt a fourth-order upwind scheme. We also use fourth-order Kreiss-Oliger dissipation to avoid high frequency noise appearing near the outer boundary. Again, we stress that in the approach adopted in this work, only the scalar field needs to be evolved numerically. The spacetime is updated following the algebraic equations from Section II.2.

Our simulations employ a logarithmic radial grid, as described in Sanchis-Gual et al. (2015b). We set the finest radial resolution close the origin, and a grid spacing of . In the first set of simulations, the outer boundary of the computational domain is placed at , far enough as to not affect the dynamics of the scalar field in the inner region during the entire simulation. The time step is chosen to which guarantees long-term stable simulations. The final time of the numerical evolutions is . In the second set of simulations, corresponding to the mirrored states, the radial extension of the numerical domain can be significantly reduced. Specifically we place the outer boundary at which allows us to use ten times less grid points than in the first setup. Finally, the time in the second set of simulations extends up to .

## Iii Results

### iii.1 Time evolution of the scalar field amplitude

#### iii.1.1 Quasi-bound frequencies

Quasi-bound states of scalar fields around BHs in the test-field approximation are configurations that have well defined (complex) frequencies. The real part gives the oscillation of the field and the imaginary part gives the rate of decay of the configuration. As mentioned in the Introduction, there is a growing body of work which has shown that for some choices of parameters quasi-bound states may survive for cosmological timescales around SMBHs and, thus, they are consistent with dark matter models based on (ultra-light) scalar fields Barranco et al. (2011, 2012, 2014); Sanchis-Gual et al. (2015a, b). Such long-lasting, quasi-bound states have also been found to exist in the nonlinear regime Okawa et al. (2014).

We turn next to discuss the effect of accretion on the quasi-bound states. We solve the Klein-Gordon system (9) using the initial data given by Eqs. (19)-(21) in the rapidly changing gravitational field of an accreting BH. As done in Sanchis-Gual et al. (2015a, b) we analyze the results of the simulations by extracting a time series for the scalar field amplitude at a set of observation points located at fixed radii (typically at ). We identify the frequencies at which the scalar field oscillates by performing a Fast Fourier transform. It is worth mentioning that the values of the frequencies do not depend on the location of the observer.

Our main results are summarized in Table 1. The first three columns report our different choices for the growth rate of the mass of the BH, , the five different scalar field masses considered in this work, and the initial amplitude of the Gaussian pulse, . The reason we vary the initial amplitude of the pulse for the different values of is to keep the initial scalar field energy (almost) constant for the different accretion rates considered, namely . The oscillation frequencies of the scalar field are reported in columns 4 and 5 for the fundamental mode and the first overtone, respectively. The imaginary part of the frequencies, shown in columns 6 to 8, is the decay rate of the oscillations of the scalar field. From left to right these three columns report the linear, quadratic, and cubic parts of analytic fits to numerical data for the different models. The last column of Table 1 also reports the mass of the BH at the end of the simulation, , computed from Eq. (18), whose evolution is displayed in Fig. 1.

Figure 2 shows the time evolution of the scalar field amplitude and spectra extracted at . In the left panels we plot the evolution for two representative scalar field masses, namely (top) and (bottom), and for the different BH mass growth rates. The criterion for the color of the different curves and its relation to follows from Fig. 1. For each panel, the decay of the amplitude of the oscillations of the scalar field is seen to depend strongly on , regardless of . For the higher growth rates (orange and blue curves), the field amplitudes by the end of the simulations have almost practically vanished. Large field amplitudes and long-lasting oscillations are most noticeable for the three lower accretion rate models of our sample. The differences for with respect to the non-accreting case become barely visible only for late times. The right panels of Fig. 2 display the corresponding power spectra. They show well-defined oscillation frequencies, best seen in the insets, of the quasi-stationary states for the models with smaller BH mass growth rate (black, red and green curves), as expected from the bandwidth theorem, while the other two models attain wider peaks.

As expected, in the non-accreting BH case, the oscillation frequencies shown in Fig. 2 and reported in Table 1 are in very good agreement with those found in Sanchis-Gual et al. (2015a) in the test-field regime and for the same scalar field masses. Despite in the current work we are not evolving the spacetime, while we did so in Sanchis-Gual et al. (2015a), the discarded back-reaction on the evolution of the scalar field is negligible because of the small energy of the scalar field in our models.

Our results also indicate that the frequencies of the quasi-bound states decrease as the mass of the BH increases. One may infer such behavior from a classical mechanics analogy. Let us consider the solution of the one-dimensional wave equation in flat spacetime with moving boundaries. If the boundary moves with a constant velocity an explicit solution can be obtained. Let us consider the oscillation of a string along the direction. The string is constrained at two points, and , and the amplitude of the wave is given by . We are interested in the solution of the equation

(22) |

in the interval , with the boundary conditions, and initial data, .

The general solution of Eq. (22) is a periodic function in with period of the form

(23) |

where and and the coefficients are determined by the initial conditions. What is important to notice in this solution is that the frequency decreases as the velocity of the boundary increases. One may think that a similar process happens when the mass of the BH increases since that implies the displacement of the position of the horizon and, thus, the displacement of the position of the boundary that confines the scalar field.

#### iii.1.2 Scalar field energy

The time evolution of the energy of the scalar field, computed with Eq. (11), is shown in solid lines in Figs. 3 and 4 for all scalar field masses, . The dashed lines in the three panels of Fig. 3 correspond to analytic fits of the numerical data, that we discuss below. These two figures show that the initial energy of the scalar field is (almost) the same for all models considered. Both figures show that the decrease of the initial energy is more significant and faster the larger the scalar field mass and the larger the BH growth.

Let us first discuss the lightest scalar field models. In figure 3, corresponding to , we see that for the smallest, non-zero value of , the decay of the scalar field energy is very close to the non-accreting case (compare red and black curves). For the higher values of the mass flux, the energy at early times begins decaying exponentially with a small slope but, at some point, the falling of scalar field energy onto the BH speeds up. The appearance of this effect depends on the BH mass growth rate, appearing earlier for higher . After some time, the process stops and the scalar field energy settles to an almost constant value, given by the energy of the part of the scalar field that has escaped to infinity. There is, however, part of the scalar field still localized around the BH in the form of quasi-stationary bound states.

On the other hand, Fig. 4 shows that for and , the decay rate of the energy of the scalar field is much faster than for the lightest models and the energy decreases almost identically regardless of the BH mass growth rate. This becomes clear if we check the imaginary part of the frequency reported in column 6 of Table 1 for these scalar field masses. The imaginary part of the frequency does not vary significantly when the growth of the BH mass is faster.

However, as we mentioned before for the smaller values of the scalar field mass displayed in Fig. 3, the energy decay changes from an exponential decay with a linear exponent (the usual feature of a quasi-stationary state around a non-accreting Schwarzschild BH) to an exponential decay with both quadratic and cubic exponents, as shown by the analytic fits. These parts are entirely due to the growth of the BH and are orders of magnitude smaller than the linear part of the exponential decay. Therefore, they are only significant for sufficiently long times. Columns 6, 7, and 8 of Table 1 report the values of the rate of decay according to the fit , where , and are the linear, quadratic and cubic coefficients, respectively. The dashed black lines in Fig. 3 fit part of the time evolution of the scalar field energy for the three models with higher BH growth rate. The values of the cubic, quadratic and linear coefficients increase with the scalar field mass, as can be inferred from Table 1.

1 | 2 | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

0.00 | 0.05 | 0.05168 | 0.05781 | 1.593E-05 | 0 | 0 | 0 | 0 | 0 | 1.000 | 8.0E05 |

5E-09 | 0.05 | 0.05168 | 0.05781 | 1.591E-05 | 9.738E-14 | … | … | … | … | 1.004 | 8.0E05 |

5E-08 | 0.05 | 0.05167 | 0.05780 | 1.590E-05 | 1.041E-12 | … | … | … | … | 1.040 | 8.0E05 |

5E-07 | 0.05 | 0.05163 | 0.05776 | 1.647E-05 | 7.612E-12 | 9.051E-18 | … | … | … | 1.221 | 4.0E05 |

5E-06 | 0.05 | 0.05137 | 0.05754 | 5.837E-06 | 6.122E-10 | -5.102E-15 | 2.179E-20 | … | … | 3.201 | 2.4E05 |

5E-05 | 0.05 | 0.05076 | 0.05688 | -2.713E-05 | 2.762E-08 | -1.928E-12 | 4.847E-17 | … | … | 6.686 | 3.8E04 |

1E-04 | 0.05 | 0.05055 | 0.05655 | 8.341E-05 | 4.054E-08 | -6.348E-12 | 3.415E-16 | … | … | 9.025 | 2.2E04 |

0.00 | 0.08 | 0.08039 | 0.08431 | 3.265E-05 | 0 | 0 | 0 | 0 | 0 | 1.000 | 4.0E05 |

5E-09 | 0.08 | 0.08039 | 0.08431 | 3.260E-05 | 2.455E-13 | … | … | … | … | 1.002 | 4.0E05 |

5E-08 | 0.08 | 0.08038 | 0.08430 | 3.272E-05 | 2.358E-12 | … | … | … | … | 1.020 | 4.0E05 |

5E-07 | 0.08 | 0.08034 | 0.08426 | 3.368E-05 | 2.042E-11 | 3.407E-17 | … | … | … | 1.221 | 4.0E05 |

5E-06 | 0.08 | 0.08012 | 0.08406 | 5.378E-05 | 2.742E-10 | -3.517E-15 | 5.500E-20 | … | … | 2.226 | 1.6E05 |

5E-05 | 0.08 | 0.07936 | 0.08337 | 2.912E-04 | -2.237E-08 | 5.539E-13 | 5.264E-17 | … | … | 3.857 | 2.7E04 |

1E-04 | 0.08 | 0.07981 | 0.08293 | -6.087E-04 | 5.364E-07 | -1.377E-10 | 1.577E-14 | -7.606E-19 | 1.332E-23 | 6.050 | 1.8E04 |

0.00 | 0.10 | 0.09975 | 0.10287 | 6.355E-05 | 0 | 0 | 0 | 0 | 0 | 1.000 | 4.0E05 |

5E-09 | 0.10 | 0.09975 | 0.10287 | 6.348E-05 | 6.340E-13 | … | … | … | … | 1.002 | 4.0E05 |

5E-08 | 0.10 | 0.09974 | 0.10286 | 6.361E-05 | 6.300E-12 | … | … | … | … | 1.020 | 4.0E05 |

5E-07 | 0.10 | 0.09969 | 0.10281 | 6.898E-05 | 2.991E-11 | 1.283E-16 | … | … | … | 1.191 | 3.5E05 |

5E-06 | 0.10 | 0.09945 | 0.10262 | 9.451E-05 | 1.602E-10 | 4.284E-15 | 9.373E-20 | … | … | 1.733 | 1.1E05 |

5E-05 | 0.10 | 0.09864 | 0.10179 | 6.781E-04 | -1.537E-07 | 1.690E-11 | -5.78E-16 | 7.136E-21 | … | 3.490 | 2.5E04 |

1E-04 | 0.10 | 0.09828 | 0.10151 | 7.818E-04 | -3.295E-07 | 6.964E-11 | -4.692E-15 | 1.153E-19 | … | 5.474 | 1.7E04 |

### iii.2 Mirrored states

In order to further understand the influence of the growth of the BH mass in the evolution of the scalar field, we turn now to describe the results corresponding to the evolution of an additional set of 21 models. These models correspond to Schwarzschild-like BH spacetimes with 7 different BH mass growth rates, namely and we impose reflecting boundary conditions for the scalar field at some radius, as done in Sanchis-Gual et al. (2016). At the location of the mirror, , and beyond, the scalar field is required to vanish. This leads to a discontinuity in the derivatives. In our simulations the mirror is located at and the boundary conditions for the scalar field are

(24) |

We consider again a Gaussian distribution for the scalar field, with and . As mentioned before, using a mirror is an idealized setup, yet it has pragmatic advantages since it allows us to perform longer runs and to study smaller growth rates. The spacetime is still given by the analytic solution of Section II.2 (hence, there are no reflections of the spacetime variables coming in from the outer boundary) and since the scalar field is now enclosed in a cavity no part of the scalar field propagates away from the BH. Therefore, all the oscillatory modes will be trapped in this case, not only those with , and the energy evolution will not be dominated by the asymptotic, scalar field energy minimum. In this idealized situation, the entire scalar field will fall into the BH with time.

The results of this new set of simulations are similar to those discussed in the preceding section. They are summarized in Table 2 for the lightest scalar field models analyzed, namely . Again, the scalar field oscillates with a fundamental frequency and higher overtones (of which Table 2 only reports the first one), but since the higher frequencies decay faster we end up with only the fundamental mode of oscillation. As for the cases without mirror, the real part of the frequency decreases with . In the case of a non-accreting BH, the decay of the energy (and of the amplitude of the oscillations) is exponential with a single linear exponent as expected. As done before, in the case of accreting BHs and long evolution times, we fit the energy decay to the analytic formula , with linear, quadratic, cubic and quartic terms. These parameters seem to depend on the BH mass growth rate and, probably, we could fit it to even higher order polynomials if we could reach even longer simulations. Indeed, some of the models reported in Table 2 can already be fitted with higher order terms (even up to 5th or 6th, coefficients and respectively) for the current evolution times considered. The 6th order term is considered for only one particular model, even if it is very small, in order to keep the highest order term always positive and, therefore, drive the fit to tend to zero for . The coefficients also depend on the scalar field mass, decreasing for smaller values of .

In Fig. 5 we plot (solid lines) time evolutions of the scalar field energy for the mirrored states of models with masses (from top to bottom) together with their corresponding fits (black dashed lines). The results discussed in the previous section become more clear in this idealized setup. Increasing the BH growth rate speeds up the decay of the energy. Moreover, the decay is longer than for the case of quasi-bound states without mirror and, contrary to what happens when there is an outgoing boundary, the energy does not relax to the value corresponding to that part of the scalar field that escapes to infinity.

### iii.3 Comparison with realistic BH mass growth rates

We can estimate the rate of the BH mass growth due to accretion using the bolometric luminosity of an Active Galaxy Nuclei (AGN) as Weihao and Yongheng (2003); Tsai et al. (2015)

(25) |

where is the bolometric luminosity, is the BH mass growth rate, is the speed of light, and is the radiative efficiency, for which we take the commonly adopted empirical value of 0.1 Yu and Tremaine (2002). Recently, bolometric luminosities larger than have been discovered Tsai et al. (2015). According to Eq. (25), this corresponds to in our units.

In figure 6 we plot several BH growth rates in geometrized units of AGN bolometric luminosities as a function of the redshift . We use data from two different samples, namely that of Woo and Urry (2002), for redshifts between 0.01 and 2.224, and the luminosities of WISE-selected galaxies reported recently in Tsai et al. (2015), which include redshifts up to . The larger luminosities are obtained precisely for the higher redshifts, corresponding to the growth phase of SMBHs. Fig. 6 shows that the maximum growth rate is only about two orders of magnitude smaller than the smallest we can afford in our simulations.

The influence of the accretion rate on the energy decay depends on the scalar field mass and decreases for smaller . Therefore, for the expected value of scalar field dark matter models, the actual value of is expected to be fairly small. However, the growth of SMBH seeds to their actual masses, , is believed to take place only during 1% of their lifetime, i.e. about years. Therefore, the role played by high order terms in the evolution of the energy could become significant during this phase and increase the amount of scalar field that may end up being absorbed by the BH during this period.

Finally, we can extrapolate our results to the realistic case of ultra-light scalar field masses. In the case of a test-field the decay rate of the dynamical resonances is related to the imaginary part of the quasiresonant frequencies, and thus its half-life time is inversely proportional to Barranco et al. (2012). For an accreting BH and the long evolution times reported in this work we have shown that terms higher than linear are important to capture the decay rate of the energy of the scalar field. We can then use the relation to solve for which corresponds to the time when the energy of the scalar field has decreased to half its initial value. The result of this exercise is depicted in Fig. 7, which shows as a function of for all five values of used in the “unmirrored” setup simulations. (We use in this figure the same color criterion for as defined in Fig. 1.) The symbols in the figure indicate the values of for the lighter values of the scalar field mass ( and 0.1) used in our simulations. Following Barranco et al. (2012) we do a least squares fit to these data and extrapolate our results to lower values of . Fig. 7 shows that even for high mass accretion rates, e.g. (red curve) quasi-bound states around accreting BHs can survive for timescales significantly longer than the age of the Universe.

## Iv Conclusions

We have studied the properties of scalar field quasi-bound states in the background of an accreting, spherically symmetric BH. We have assumed that the BH mass grows due to matter accretion, describing the effect of the increase of the BH mass on the properties of the quasi-bound states of the surrounding scalar field distribution. For our study we have solved the Klein-Gordon equation numerically, mimicking the evolution of the spacetime through a sequence of exact Schwarzschild BH solutions of increasing mass, using analytic expressions which depend only on the BH mass parameter. Our numerical approach has been limited to spherical symmetry and has been based on spherical coordinates and a PIRK numerical scheme for the time update of the Klein-Gordon equation. To study the effect of accretion on the scalar field evolution in affordable computational times, we have resorted to large growth rates, at best two orders of magnitude larger than the largest observational estimations.

By performing a Fourier transform of the time series of our numerical data we have been able to characterize the scalar field states by their distinctive oscillation frequencies. It has been found that the frequencies decrease with increasing BH mass. Moreover, accretion results in a significative increase of the exponential decay of the scalar field energy due to the presence of terms of order higher than linear in the exponent. These terms are zero in the non-accreting Schwarzschild BH case, resulting in a linear exponential fit. Our computational setup has considered both outgoing and reflecting boundary conditions at the outer radial boundary, the latter describing a scalar field enclosed in a cavity. By imposing reflecting boundary conditions at a finite distance the scalar field does not escape to infinity and we can isolate the influence of accretion. Such configuration, albeit artificial, has helped us to study the higher order terms that appear in the exponential decay of the energy. Finally, we have compared our BH mass growth rates with estimates from observational surveys, and we have been able to extrapolate our results to realistic values of the scalar field mass . We have found that even for the high mass accretion rates considered in this work, e.g. , quasi-bound states around accreting BHs can survive for cosmological timescales. The results obtained in this paper add further support to the intriguing possibility of the existence of dark matter halos based on (ultra-light) scalar fields surrounding SMBHs at galactic centers.

###### Acknowledgements.

This work was supported in part by the Spanish MINECO (AYA2013-40979-P), by the Generalitat Valenciana (PROMETEOII-2014-069), by the CONACyT-México, ICF-UNAM, and by the Max-Planck-Institut für Astrophysik. The computations have been performed at the Servei d’Informàtica de la Universitat de València.## References

- Ferrarese and Ford (2005) L. Ferrarese and H. Ford, Space Science Reviews 116, 523 (2005), eprint astro-ph/0411247.
- Kormendy and Ho (2013) J. Kormendy and L. C. Ho, Annual Review Astron. Astrophys. 51, 511 (2013), eprint 1304.7762.
- Kauffmann and Haehnelt (2000) G. Kauffmann and M. Haehnelt, Monthly Not. R. Astron Soc. 311, 576 (2000), eprint astro-ph/9906493.
- Volonteri et al. (2003) M. Volonteri, F. Haardt, and P. Madau, Astrophys.J. 582, 559 (2003), eprint astro-ph/0207276.
- Barth et al. (2003) A. J. Barth, P. Martini, C. H. Nelson, and L. C. Ho, Astrophys. J. 594, L95 (2003), eprint astro-ph/0308005.
- Willott et al. (2005) C. J. Willott, W. J. Percival, R. J. McLure, D. Crampton, J. B. Hutchings, M. J. Jarvis, M. Sawicki, and L. Simard, Astrophys.J. 626, 657 (2005), eprint astro-ph/0503202.
- Mortlock et al. (2011) D. J. Mortlock, S. J. Warren, B. P. Venemans, M. Patel, P. C. Hewett, R. G. McMahon, C. Simpson, T. Theuns, E. A. Gonzáles-Solares, A. Adamson, et al., Nature 474, 616 (2011).
- Volonteri (2010) M. Volonteri, The Astronomy and Astrophysics Review 18, 279 (2010).
- Madau and Rees (2001) P. Madau and M. J. Rees, The Astrophysical Journal Letters 551, L27 (2001).
- Barausse et al. (2015) E. Barausse, J. Bellovary, E. Berti, K. Holley-Bockelmann, B. Farris, B. Sathyaprakash, and A. Sesana, Journal of Physics Conference Series 610, 012001 (2015), eprint 1410.2907.
- Alexander and Natarajan (2014) T. Alexander and P. Natarajan, Science 345, 1330 (2014).
- Bromm and Loeb (2003) V. Bromm and A. Loeb, The Astrophysical Journal 596, 34 (2003).
- Haiman and Loeb (2001) Z. Haiman and A. Loeb, The Astrophysical Journal 552, 459 (2001).
- Begelman et al. (2006) M. C. Begelman, M. Volonteri, and M. J. Rees, Monthly Notices of the Royal Astronomical Society 370, 289 (2006).
- Begelman (2010) M. C. Begelman, Monthly Notices of the Royal Astronomical Society 402, 673 (2010).
- Mayer et al. (2010) L. Mayer, S. Kazantzidis, A. Escala, and S. Callegari, Nature 466, 1082 (2010).
- Barranco et al. (2011) J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, M. Megevand, M. Alcubierre, D. Núñez, and O. Sarbach, Phys.Rev.D 84, 083008 (2011), eprint 1108.0931.
- Barranco et al. (2012) J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, M. Megevand, et al., Phys.Rev.Lett. 109, 081102 (2012), eprint 1207.2153.
- Barranco et al. (2014) J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, M. Megevand, et al., Phys.Rev. D89, 083006 (2014), eprint 1312.5808.
- Sanchis-Gual et al. (2015a) N. Sanchis-Gual, J. C. Degollado, P. J. Montero, and J. A. Font, Phys. Rev. D 91, 043005 (2015a), eprint 1412.8304.
- Sanchis-Gual et al. (2015b) N. Sanchis-Gual, J. C. Degollado, P. J. Montero, J. A. Font, and V. Mewes, Phys. Rev. D 92, 083001 (2015b), eprint 1507.08437.
- Hannam et al. (2007) M. Hannam, S. Husa, D. Pollney, B. Brügmann, and N. O. Murchadha, Phys. Rev. Lett. 99, 241102 (2007), URL http://link.aps.org/doi/10.1103/PhysRevLett.99.241102.
- Baumgarte and Naculich (2007) T. W. Baumgarte and S. G. Naculich, Phys. Rev. D 75, 067502 (2007), URL http://link.aps.org/doi/10.1103/PhysRevD.75.067502.
- Willott et al. (2010) C. J. Willott, L. Albert, D. Arzoumanian, J. Bergeron, D. Crampton, P. Delorme, J. B. Hutchings, A. Omont, C. Reylé, and D. Schade, The Astronomical Journal 140, 546 (2010).
- Fan et al. (2006) X. Fan, M. A. Strauss, R. H. Becker, R. L. White, J. E. Gunn, G. R. Knapp, G. T. Richards, D. P. Schneider, J. Brinkmann, and M. Fukugita, The Astronomical Journal 132, 117 (2006).
- Springel et al. (2005) V. Springel, T. Di Matteo, and L. Hernquist, Monthly Notices of the Royal Astronomical Society 361, 776 (2005).
- Hu et al. (2006) J. Hu, Y. Shen, Y.-Q. Lou, and S. Zhang, Monthly Notices of the Royal Astronomical Society 365, 345 (2006).
- Di Matteo et al. (2008) T. Di Matteo, J. Colberg, V. Springel, L. Hernquist, and D. Sijacki, The Astrophysical Journal 676, 33 (2008).
- Tsai et al. (2015) C.-W. Tsai, P. R. Eisenhardt, J. Wu, D. Stern, R. J. Assef, A. W. Blain, C. R. Bridge, D. J. Benford, R. M. Cutri, R. L. Griffith, et al., The Astrophysical Journal 805, 90 (2015).
- Sanchis-Gual et al. (2016) N. Sanchis-Gual, J. C. Degollado, P. J. Montero, J. A. Font, and C. Herdeiro, Physical Review Letters 116, 141101 (2016).
- Montero and Cordero-Carrion (2012) P. J. Montero and I. Cordero-Carrion, Phys.Rev. D85, 124037 (2012), eprint 1204.5377.
- Cordero-Carrión and Cerdá-Durán (2012) I. Cordero-Carrión and P. Cerdá-Durán, ArXiv e-prints (2012), eprint 1211.5930.
- Cordero-Carrión and Cerdá-Durán (2014) I. Cordero-Carrión and P. Cerdá-Durán, Advances in Differential Equations and Applications, SEMA SIMAI Springer Series Vol. 4 (Springer International Publishing Switzerland, Switzerland, 2014).
- Okawa et al. (2014) H. Okawa, H. Witek, and V. Cardoso, Phys.Rev. D89, 104032 (2014), eprint 1401.1548.
- Weihao and Yongheng (2003) B. Weihao and Z. Yongheng, arXiv preprint astro-ph/0305095 (2003).
- Yu and Tremaine (2002) Q. Yu and S. Tremaine, Monthly Notices of the Royal Astronomical Society 335, 965 (2002).
- Woo and Urry (2002) J.-H. Woo and C. M. Urry, The Astrophysical Journal 579, 530 (2002).