# Primordial black holes as dark matter: constraints from compact ultra-faint dwarfs

Primordial black holes as dark matter: constraints from compact ultra-faint dwarfs Abstract The ground-breaking detections of gravitational waves from black hole mergers by LIGO have rekindled interest in primordial black holes (PBHs) and the possibility of dark matter being composed of PBHs. It has been suggested that PBHs of tens of solar masses could serve as dark matter candidates. Recent analytical studies demonstrated that compact ultra-faint dwarf galaxies can serve as a sensitive test for the PBH dark matter hypothesis, since stars in such a halo-dominated system would be heated by the more massive PBHs, their present-day distribution can provide strong constraints on PBH mass. In this study, we further explore this scenario with more detailed calculations, using a combination of dynamical simulations and Bayesian inference methods. The joint evolution of stars and PBH dark matter is followed with a Fokker–Planck code phaseflow. We run a large suite of such simulations for different dark matter parameters, then use a Markov chain Monte Carlo approach to constrain the PBH properties with observations of ultra-faint galaxies. We find that two-body relaxation between the stars and PBH drives up the stellar core size, and increases the central stellar velocity dispersion. Using the observed half-light radius and velocity dispersion of stars in the compact ultra-faint dwarf galaxies as joint constraints, we infer that these dwarfs may have a cored dark matter halo with the central density in the range of 1–2 M⊙pc − 3, and that the PBHs may have a mass range of 2–14 M⊙ if they constitute all or a substantial fraction of the dark matter. methods: numerical, dark matter 1 INTRODUCTION In the standard lambda cold dark matter (ΛCDM) cosmology, the nature of the dark matter (DM) remains elusive. The possibility of DM being composed of primordial black holes (PBHs; Hawking 1971), arising naturally from certain inflation models (e.g. Frampton et al. 2010), has been investigated in the context of galaxy formation early on (e.g. Carr & Hawking 1974; Meszaros 1975), and it has not yet been completely ruled out (see recent reviews by Khlopov 2010; Carr, Kühnel & Sandstad 2016). Recently, it was suggested that three mass windows, around 5 × 10 − 16 M⊙, 2 × 10 − 14 M⊙, and 25–100 M⊙, respectively, are still open for DM consists of PBHs (Carr et al. 2017). The recent detections of gravitational waves by LIGO are believed to originate from merging black holes with masses of order 10–50 M⊙, significantly higher than known BHs in Galactic X-ray binaries (The LIGO Scientific Collaboration et al. 2016a,b, 2017a,b). These discoveries have renewed theoretical interest in PBHs as a DM candidate. It was suggested by Bird et al. (2016) that PBHs with a mass around 30 M⊙ are loosely constrained by micro-lensing experiments and that the merger rate of binary PBHs is consistent with the LIGO detections. If true, PBHs in this mass range could be an attractive DM candidate to solve some well-known problems on the galactic scale (see Weinberg et al. 2015), in particular the ‘cusp versus core problem’ (hereafter, we refer to dark matter made of primordial black holes as PBH-DM). The PBHs would produce DM cores instead of cusps in the central galaxy because any ‘temperature inversion’ present in a cuspy profile will lead to a removal of the density cusp by increasing the velocity dispersion (Quinlan 1996). Solutions to the ‘cusp versus core problem’ based on baryonic heating of DM are shown to be mass-dependent (Di Cintio et al. 2014; Pontzen & Governato 2014; Chan et al. 2015). The survival of a globular cluster in the ultra-faint dwarf galaxy Eridanus II, however, hints to a DM core in a galaxy halo where any baryonic effect would be weak (Amorisco 2017; Contenta et al. 2017). In addition, a cosmology model with PBHs as DM indicates a rather different evolution picture from that of other types of DM, such as warm dark matter, fuzzy dark matter (Hu, Barkana & Gruzinov 2000; Hui et al. 2017), and self-interacting dark matter (Vogelsberger, Zavala & Loeb 2012). Contrary to a common feature of suppressed small-scale perturbations in these models (Vogelsberger et al. 2016), fluctuations in the number density of PBH impose additional small-scale power (Afshordi, McDonald & Spergel 2003). As a result, low-mass PBH-DM haloes will form ahead of ΛCDM haloes, while in warm/fuzzy dark matter, the DM halo formation is considerably delayed, which has a strong consequence on the onset of reionization (Yoshida et al. 2003; Hirano, Sullivan & Bromm 2017). In addition to the increased collapsed halo at high redshift, Kashlinsky (2016) has shown that gas accretion on to the PBHs could explain the observed high coherence between cosmic infrared and the soft-X-ray background. Despite these appealing features, it remains unclear what fraction of the DM can PBHs constitute, and what mass range can these PBHs have. Fortunately, a population of compact (with projected half-light radius of ∼30 pc), ultra-faint (∼1000 L⊙) dwarf galaxies has recently been discovered (e.g. Koposov et al. 2015; Bechtol et al. 2015), and they are expected to provide stringent constraints on the PBH-DM. Using analytical approaches, Brandt (2016) suggested that two-body relaxation between stars and PBH-DM would drive a substantial size evolution in galaxies, which is inconsistent with the observed small size of the ultra-faint dwarfs, while Koushiappas & Loeb (2017) argued that mass segregation would rearrange the stars into an otherwise detectable ring in the surface brightness profile of Segue I. In order to improve these analytical studies, we combine two methods, Fokker–Planck (FP) simulations and Markov chain Monte Carlo (MCMC) modelling, to study the evolution of galaxy haloes containing both stars and PBH-DM and to compare with observations. We first use an FP code to perform a large set of galaxy simulations with different parameters for the PBH-DM to follow the interaction and evolution of the stars and PBHs, then we employ an MCMC method to constrain the PBH properties using the observed half-light radius and velocity dispersion of stars in compact ultra-faint galaxies. Our modelling presents a significant improvement over the previous studies based on an analytical form of energy diffusion coefficient and other simplifying assumptions (Brandt 2016; Koushiappas & Loeb 2017). This paper is organized as follows. In Section 2, we describe the numerical code and MCMC algorithm used in our simulations and analysis. In Section 3, we present the dynamical evolution of the stars and PBHs in the halo systems, and constraints from observations. We discuss the implications of our study and its limitations in Section 4, and summarize our main findings in Section 5. 2 METHODS 2.1 Numerical code for dynamical simulations In this study, we use an FP code, phaseflow1 (Vasiliev 2017) to follow the dynamical evolution of stars and PBH-DM in the galaxy haloes. phaseflow solves the one-dimensional FP equation for the distribution function in energy, using a high-accuracy finite element method. The code is well tested with problems such as core collapse and the formation of Bahcall–Wolf cusp around a central massive black hole. phaseflow can handle multiple mass components, which enables our calculations of the evolution of a halo of stars and DM consisting of PBHs. Once the initial conditions of the two components are specified, phaseflow evolves the distribution function of all components while at the same time updating the gravitational potential. System diagnosis such as density profile and velocity dispersions are computed automatically, which greatly simplifies our analysis. Throughout the paper, we set the Coulomb logarithm (Binney & Tremaine 2008) ln Λ to be 15. The stars in the DM halo are modelled with a Plummer sphere (Plummer 1911):   \begin{eqnarray} \rho (r) = \frac{3M_{*}}{4\pi R_{0,*}^3} \bigg (1+\frac{r^2}{R_{0,*}^2} \bigg )^{-5/2}, \end{eqnarray} (1)where the stellar system is characterized by the total stellar mass M* and its scale radius R0, *, respectively. In this calculation, we set M* = 103 M⊙, and each star particle has an equal mass of 1 M⊙. The extended PBH-DM halo follows a Dehnen sphere (Dehnen 1993; Contenta et al. 2017):   \begin{eqnarray} \rho (r) = \frac{(3-\gamma ){M_{\rm DM}}}{4\pi R_{\rm 0, DM}^3} \bigg (\frac{r}{R_{\rm 0, DM}}\bigg )^{-\gamma } \bigg (1 + \frac{r}{R_{\rm 0, DM}}\bigg )^{\gamma -4}, \end{eqnarray} (2)where MDM is the total mass and R0, DM the scale radius of the halo, respectively. The Dehnen sphere is a special case of the more general α–β–γ density profile (Zhao 1996), where the density distribution is determined by the inner slope γ and the outer slope β and the steepness parameter α, which control how quickly the inner slope transits into the outer slope. Note that when γ = 1, the density distribution becomes the Hernquist profile (Hernquist 1990). In order to compare with previous studies, we also consider the Navarro–Frenk–White (NFW) profile (Navarro, Frenk & White 1997) in our simulations. 2.2 Markov chain Monte Carlo for parameter constraints 2.2.1 Observational data and the likelihood function The MCMC is an effective method to explore multiple key parameters in a complex system. In our model, we use both the stellar half-light radius and the stellar velocity dispersion of the compact ultra-faint dwarfs to constrain the parameter space, which requires ∼105 MCMC realizations. Fortunately, phaseflow is an extremely efficient code, and it takes only several seconds to run a single model. Thus, we combine a large grid of FP simulations performed using phaseflow and MCMC calculations using the emcee code. The emcee code uses an affine-invariant ensemble sampler to efficiently random walk the parameter space (Foreman-Mackey et al. 2013). The likelihood function of our modelling is described in the following form:   \begin{eqnarray} \mathcal {L} = \prod _{i=1}^{N} \mathcal {N} (r_{\rm h, M}, r_{\rm h, O}, \epsilon _{\rm h, O})\, \mathcal {N} (\sigma _{\rm M}, \sigma _{\rm O}, \epsilon _{\rm \sigma , O}), \end{eqnarray} (3)where $$\mathcal {N}$$ is a normal distribution, N is the sample size, rh, M is the model output of 3D half-mass radius, rh, O, and σO are the observed half-light radius and velocity dispersion, respectively. The associated uncertainties are εh, O and εσ, O. We use the top five ultra-faint dwarf galaxies compiled by Brandt (2016), which have measurements of both size and stellar velocity dispersion available. The data are repeated here in Table 1 for convenience but we refer the readers to Brandt (2016) for details and references of each measurement. Note that there are discrepancies in the reported half-mass radius of the two newly discovered Ret II and Hor I between Koposov et al. (2015) and Bechtol et al. (2015). In each case, we have chosen the smaller value of the two measurements. The reason for such a sample is based on the implicit assumption (inherited from Brandt 2016) that these five galaxies form a distinct class of faint dwarfs with similar size and velocity dispersion, which can be compared to a single theoretical model. Lastly, we apply a factor of 1.3 to de-project the observed half-mass radius to 3D following Wolf et al. (2010). Table 1. A sample of observed compact ultra-faint dwarf galaxies used in our modelling from Brandt (2016). Galaxy name  Projected rh (pc)  σ* (km s−1)  LV (L⊙)  Wil I  25 ± 6  $$4.3^{+2.3}_{-1.3}$$  1000  Seg I  $$29^{+8}_{-5}$$  $$3.9^{+0.8}_{-0.8}$$  300  Seg II  35 ± 3  $$3.4^{+2.5}_{-1.2}$$  900  Ret II  $$32^{+2}_{-1}$$  $$3.2^{+1.6}_{-0.5}$$  1500  Hor I  $$25^{+9}_{-4}$$  $$4.9^{+2.8}_{-0.9}$$  2000  Galaxy name  Projected rh (pc)  σ* (km s−1)  LV (L⊙)  Wil I  25 ± 6  $$4.3^{+2.3}_{-1.3}$$  1000  Seg I  $$29^{+8}_{-5}$$  $$3.9^{+0.8}_{-0.8}$$  300  Seg II  35 ± 3  $$3.4^{+2.5}_{-1.2}$$  900  Ret II  $$32^{+2}_{-1}$$  $$3.2^{+1.6}_{-0.5}$$  1500  Hor I  $$25^{+9}_{-4}$$  $$4.9^{+2.8}_{-0.9}$$  2000  View Large 2.2.2 Selections of prior distributions Total dark matter mass MDM: Compact ultra-faint dwarf galaxies from our sample have the least amount of stars (∼1000 L⊙) ever known, but their total mass remains unknown. It was suggested that they may reside in the least massive haloes where atomic cooling is efficient, and that their total mass is around 109 M⊙ (Sawala et al. 2015; Wheeler et al. 2015; Zhu et al. 2016), since such a halo mass is consistent with the M*–MDM relations at the low mass end (Garrison-Kimmel et al. 2017). Recently, Ma et al. (2017) resolved the formation of such galaxies in their simulations, and reported that the least massive halo, which has a mass of 108M⊙ at z = 5, had already stopped star formation since z = 8. On the other hand, Read, Agertz & Collins (2016) modelled isolated dwarf galaxies and found that haloes with a total mass of 108M⊙ were able to match several dwarf satellites of the Milky Way. Therefore, we choose a weakly informative prior of MDM, which follows a log-normal distribution:   \begin{eqnarray} \log ({M_{\rm DM}}) \sim \mathcal {N}(9, 0.5). \end{eqnarray} (4) The choice of this prior is well justified in that the physics involved is relatively straightforward. These low-mass galaxies must be just massive enough, above the thresholds set by atomic cooling and UV background radiation, for stars to form (e.g. Okamoto, Gao & Theuns 2008; Gnedin & Kaurov 2014). The width of the distribution also ensures that haloes down to 2–3 × 108 M⊙ are covered. DM scale radius R0, DM:From the discussion below (Section 3.1), we argue that the inner density profile of the halo will form a core (γ = 0) due to relaxation effect. However, the size of the scale radius needs to be treated as another free parameter in our model. This is largely due to a lack of self-consistent collisional N-body simulations from cosmological initial conditions. The assumption that the characteristic scales of the PBH-DM haloes are similar to the ΛCDM models could lead to serious systematics. Thus, we adopt a log-uniform distribution of R*, DM for the range from 100 to 2000 pc:   \begin{eqnarray} \log (R_{\rm *, DM}) \sim \mathcal {U}(\log (100), \log (2000)). \end{eqnarray} (5) Stellar scale radius R0, *:The heating of stars due to two-body relaxation increases the scale radius of stellar distribution roughly as R* ∼ t0.5 (Brandt 2016), which indicates that the present-day size does not depend sensitively on the initial value of R0, *. We thus use a log-uniform distribution for R0, * for the range from 10 to 50 pc.   \begin{eqnarray} \log (R_{\rm 0, *}) \sim \mathcal {U}(\log (10), \log (50)). \end{eqnarray} (6) The lower limit corresponds to the largest Galactic globular clusters and the upper limit corresponds to the average size of the currently observed compact ultra-faint dwarfs. Primordial BH mass MBH: Lastly, the mass of the PBHs is not well constrained (Carr et al. 2017). PBHs above 100 M⊙ can be ruled out by CMB observation and the Galactic wide binaries (see Carr et al. 2016). However, since we are most interested in the range of [1, 100] M⊙, we assume a log-uniformly distribution for MBH in this range. 3 RESULTS 3.1 Cusp-to-core transition in PBH-DM haloes from relaxation effect In a CDM halo, the density profile usually shows a central density cusp with a slope in the range of 0 < γ < 2, and the DM velocity dispersion peaks at the scale radius. However, if the DM is composed of PBHs, two-body relaxation will soften the central cusp, and as a result of the collisional heating, the velocity dispersion in the central region will increase while the density will drop at the same time. This process leads to the rapid formation of a core (Quinlan 1996) as a result of ‘temperature inversion’, in which the colder dense cusp will be heated by a hotter envelope. Fig. 1 illustrates how this relaxation effect transforms a PBH-DM halo from a γ = 1 cusp to a γ = 0 core. The DM halo has a mass of 2 × 109 M⊙ composed of 30 M⊙ PBHs, and an initial cuspy density profile of γ = 1 and R0, DM = 500 pc. Collisional relaxation quickly removes the central density cusp and increases the central velocity dispersion. After 12 Gyr, this PBH-DM halo can be approximated by a less massive (109M⊙) one with a core density profile of γ = 0 and a characteristic scale radius of R0, DM = 160 pc. In addition, the evolution of an NFW halo with a scale radius of 300 pc and truncated at 10 kpc is shown in Fig. 1 for comparison. The scale radius is chosen to match the initial cuspy Dehnen profile. Within the central 200 pc, both the density and the velocity dispersion profiles of the NFW halo match perfectly with those of the cuspy Dehnen model. While the NFW profile leads different asymptotic behaviour towards large r, but the dynamical region we focus in this study is within the central ∼ 100 pc, where both NFW and the cuspy Dehnen profiles produce the same results. Figure 1. View largeDownload slide The transition from cusp to core (top panel) and the increase of central velocity dispersion (bottom panel) due to collisional relaxation effect in a PBH-DM halo. The DM halo has a mass of 2 × 109M⊙ consisting of 30 M⊙ PBHs. It has an initial density profile of γ = 1 (Hernquist profile) and R0, DM = 500 pc, and a diminishing velocity dispersion as r → 0, as represented by the thick black curves in both panels. The thin lines show the evolution of the density profile and the velocity dispersion from t = 0 to 12 Gyr with an interval of 1.5 Gyr. The red dashed lines shows the best matching PBH-DM halo of 109M⊙ with a core density profile of γ = 0 and R0, DM = 160 pc, and an increased central velocity dispersion. For comparison, we also show the evolution of an NFW halo, which is truncated at 10 kpc, at t = 0 and 12 Gyr, respectively (green dashed curves). Figure 1. View largeDownload slide The transition from cusp to core (top panel) and the increase of central velocity dispersion (bottom panel) due to collisional relaxation effect in a PBH-DM halo. The DM halo has a mass of 2 × 109M⊙ consisting of 30 M⊙ PBHs. It has an initial density profile of γ = 1 (Hernquist profile) and R0, DM = 500 pc, and a diminishing velocity dispersion as r → 0, as represented by the thick black curves in both panels. The thin lines show the evolution of the density profile and the velocity dispersion from t = 0 to 12 Gyr with an interval of 1.5 Gyr. The red dashed lines shows the best matching PBH-DM halo of 109M⊙ with a core density profile of γ = 0 and R0, DM = 160 pc, and an increased central velocity dispersion. For comparison, we also show the evolution of an NFW halo, which is truncated at 10 kpc, at t = 0 and 12 Gyr, respectively (green dashed curves). In order to determine the time-scale of the transition from cusp-to-core profiles, we set up cuspy DM haloes consisting of 30 M⊙ PBHs with a total mass in the range from 105 to 109 M⊙ at redshift z = 10, based on the mass-concentration relation from Diemer & Kravtsov (2015). We evolve the PBH-DM haloes using the FP code phaseflow. These simulations show that the collisional relaxation effect removes the central cusp almost instantaneously, and the DM core grows quickly in size. For instance, it takes only ∼0.05 Gyr for all the haloes to develop a core of ∼10 pc, and ∼0.18 Gyr to reach ∼20 pc. If formed via hierarchical assembly, it would be difficult for DM haloes of $$10^8\text{--}10^9\, \rm {M_{\rm {\odot }}}$$ to retain density cusps in the first place since the central density of their less massive progenitors has already been lowered. We speculate that the density profile of a PBH-DM halo, similar to self-interacting DM (Vogelsberger et al. 2012), could be described by a truncated singular sphere with a sizable core (Shapiro, Iliev & Raga 1999). In the absence of a self-consistent treatment of density profile of PBH-DM haloes, we thus use a cored density profile with γ = 0 for any initial DM scale radius R0, DM as our default choice for simulations in the following sections. In the Discussion Section, we will include a test using γ = 1 to examine the outcome with cuspy density profiles. 3.2 Stars in PBH-DM haloes In a system with two components of different masses, two-body encounters and energy exchange would lead to change in the distribution of density and energy. To investigate the evolution of a DM halo with stars and PBHs, we set up a DM halo with a mass of 2 × 109 M⊙, which contains 103 M⊙ of stars, each being 1 M⊙, and the DM is composed of 30 M⊙ PBHs. The Plummer sphere for the stars has a total mass of 103M⊙ and a scale radius of 5 pc, while the Dehnen sphere for the DM has a total mass of 2 × 109M⊙ and a scale radius R0, DM = 500 pc. The system is then followed with the FP code phaseflow, and the results are shown in Fig. 2. Due to the huge difference in both mass and size, the two components have different initialization in both density profile and velocity dispersion, and strong two-body relaxation and energy exchange leads to significant change in the stellar density profile and velocity dispersion. Because the PBH is much more massive than the star, heating from PBHs on the stars is the dominant driver for the evolution of the stellar component. After 12 Gyr, the central stellar density drops by nearly 3 orders of magnitude, as shown in the upper panel of Fig. 2, while the central stellar velocity dispersion increases by a factor of ∼8, as in the lower panel. Figure 2. View largeDownload slide Effects of two-body relaxation and energy exchange on the density profile (top panel) and velocity dispersion (bottom panel) of a two-component halo system with stars and PBH-DM. The halo has a total DM mass of 2 × 109 M⊙ composed of 30 M⊙ PBHs, and a total stellar mass of 103 M⊙ consisting of 1 M⊙ stars. The blue curves represent the PBH-DM component at t = 0 (dashed lines) and t = 12 Gyr (solid lines), while the red curves represent the stellar components at t = 0 (dashed lines) and t = 12 Gyr (solid lines). Figure 2. View largeDownload slide Effects of two-body relaxation and energy exchange on the density profile (top panel) and velocity dispersion (bottom panel) of a two-component halo system with stars and PBH-DM. The halo has a total DM mass of 2 × 109 M⊙ composed of 30 M⊙ PBHs, and a total stellar mass of 103 M⊙ consisting of 1 M⊙ stars. The blue curves represent the PBH-DM component at t = 0 (dashed lines) and t = 12 Gyr (solid lines), while the red curves represent the stellar components at t = 0 (dashed lines) and t = 12 Gyr (solid lines). As the density profile of the stars slowly diffuses out, the half-light radius increases. This process depends on the density and velocity dispersion profile of the PBH halo. We note that in Fig. 2, the stellar velocity dispersion is well below that of PBH-DM at t = 12 Gyr. Hence, equipartition of energy between the two mass species, i.e. σ(m)2 ∝ m−1, has not been achieved. As first pointed out by Spitzer (1969), full energy equipartition is possible only if the mass fraction of heavy species is below some critical value. In our setup, the total mass ratio between PBHs and stars is well above the critical value ∼0.001(= 0.16(1/30)1.5). As a result, the central regime will be devoid of stars, leading to only a ‘partial equipartition’. Recently, Trenti & van der Marel (2013) and Bianchini et al. (2016) reported partial equipartition in their direct N-body simulations of globular clusters. Of course, our example is quite extreme since the mass density profile is completely dominated by PBHs at almost all radii. The size increase due to heating of PBH appears to be slightly slower than the analytical result obtained by Brandt (2016). To further investigate this, we performed a set of FP simulations with phaseflow for different initial stellar density distributions, by varying R0, * from 3 to 15 pc, within a fixed Dehnen sphere with a total DM mass of 2 × 109M⊙ and R0, DM = 500 pc. The resulting size evolution of the different stellar components is shown in Fig. 3, in comparison with that of Brandt (2016). Our simulations show a slower growth rate, rh ∝ t0.4, than that of rh ∝ t0.5 by Brandt (2016). The heating rate of a less concentrated stellar component is also slower than those of more concentrated ones. For example, the magenta line on top of Fig. 3 shows the size evolution of a stellar core from R0, * = 15 pc, which only approaches the asymptotic t0.4 trend by the end of the integration. Our results suggest that the final size of the stellar core only depends weakly on the initial size, and that after ∼10  Gyr, a two-component PBH-DM halo of 2 × 109M⊙ would produce a stellar core of ∼50 pc, regardless of its initial size. Figure 3. View largeDownload slide The evolution of 3D half-light radius of stellar components with various initial conditions in a PBH-DM halo from our simulations, in comparison with analytical result from Brandt (2016). The different colour represents different initial condition with R0, * varying from 3, 6, 9, 12, and 15 pc, respectively. Our results show an rh ∼ t0.4 size growth rate, slower than rh ∼ t0.5 reported by Brandt (2016). Figure 3. View largeDownload slide The evolution of 3D half-light radius of stellar components with various initial conditions in a PBH-DM halo from our simulations, in comparison with analytical result from Brandt (2016). The different colour represents different initial condition with R0, * varying from 3, 6, 9, 12, and 15 pc, respectively. Our results show an rh ∼ t0.4 size growth rate, slower than rh ∼ t0.5 reported by Brandt (2016). 3.3 Parameter space of PBH-DM haloes from FP simulations In order to explore the parameter space of DM haloes of the observed compact ultra-faint dwarfs in Table 1, we use the FP code phaseflow to perform a grid of dynamical simulations of the two-component systems consisting of stars and PBH-DM with different halo mass and size. We vary the total DM mass from 108 to 5 × 109 M⊙, and vary the scale radius of the DM Dehnen sphere from 100 to 1000 pc. We consider 50 values for each of the two DM halo parameters, MDM and R0, DM, resulting in a total of 2500 simulations. We assume the DM is composed of PBHs of 30 M⊙, a total stellar mass of 103 M⊙, and an initial stellar scale radius of R0, * = 20 pc. Each system is integrated for a duration of 12 Gyr. The final half-mass radius and the projected central velocity dispersion of the stars are shown in Fig. 4, as functions of the DM parameter space of total DM masses MDM and scale radius R0, DM. The shaded regions in both the top and bottom panels indicate the allowed range of DM parameters from the observations of the compact ultra-faint dwarfs. From the stellar half-mass radius distribution in the top panel, given the same total DM mass, the spreading of the stellar component is more efficient when the scale radius is larger. While with the same scale radius, less massive haloes show more efficient heating. When the halo is less dense (either due to larger scale radius or lower mass), the initial velocity dispersions of both stars and PBHs are also much lower. In this case, the relative increase of stellar velocity dispersion is larger than in a more dense halo. The preferred range of rh ∼ 25–40 pc from observations demand a narrow range of R0, DM < 400 pc and MDM < 5 × 109 M⊙. On the other hand, from the stellar velocity dispersion in the bottom panel, the preferred range of σ ∼ 3–5 km s−1 dictates a range of 300 < R0, DM < 1000 pc and MDM > 109 M⊙. Figure 4. View largeDownload slide The distribution of 3D half-mass radius (top panel) and 1D velocity dispersion (bottom panel) of the stellar components in the PBH-DM halo parameter space, from 2500 FP simulations using the phaseflow code. There are 50 grids for each of the DM halo parameters, total DM mass MDM and scale radius R0, DM. The DM is assumed to consist of PBHs of 30 M⊙. The allowed regions by observations are highlighted and hatched in both panels which, however, have little overlap with each other. Figure 4. View largeDownload slide The distribution of 3D half-mass radius (top panel) and 1D velocity dispersion (bottom panel) of the stellar components in the PBH-DM halo parameter space, from 2500 FP simulations using the phaseflow code. There are 50 grids for each of the DM halo parameters, total DM mass MDM and scale radius R0, DM. The DM is assumed to consist of PBHs of 30 M⊙. The allowed regions by observations are highlighted and hatched in both panels which, however, have little overlap with each other. Clearly, it is difficult to reconcile the two constraints simultaneously. A smaller stellar size would require a DM halo with smaller scale radius, thus reducing the heating from relaxation by increasing the velocity dispersion. Given PBH mass of 30 M⊙ and a typical size of ∼40 pc (de-projected), there is little room in the plane of MDM and r0, DM to meet both the constraints. 3.4 MCMC constraints on PBH-DM from compact ultra-faint dwarfs In order to constrain the PBH-DM parameters from observations, we use the stellar half-mass radius and stellar velocity dispersion of the observed compact ultra-faint dwarfs in Table 1 as priors and perform a total of 105 MCMC realizations using the emcee code. The resulting posterior probability distributions of the four model parameters, the DM scale radius R0, DM, initial stellar component scale radius R0, *, total DM mass MDM, and PBH mass MBH, are shown in Fig. 5, respectively. Figure 5. View largeDownload slide Posterior probability distributions of the four PBH-DM parameters constrained by the stellar half-mass radius and stellar velocity dispersion of observed compact ultra-faint dwarfs: the DM scale radius R0, DM, initial stellar component scale radius R0, *, total DM mass MDM, and PBH mass MBH. The contours indicate the 1σ and 2σ uncertainty region. This plot is generated using the python package ChainConsumer (Hinton 2016). Figure 5. View largeDownload slide Posterior probability distributions of the four PBH-DM parameters constrained by the stellar half-mass radius and stellar velocity dispersion of observed compact ultra-faint dwarfs: the DM scale radius R0, DM, initial stellar component scale radius R0, *, total DM mass MDM, and PBH mass MBH. The contours indicate the 1σ and 2σ uncertainty region. This plot is generated using the python package ChainConsumer (Hinton 2016). The posterior surface in the plane of R0, DM and M0, DM can be understood with the help of Fig. 4, as a balance between the requirements of small sizes and small velocity dispersions. The contours in the plane of R0, DM and M0, DM roughly cover a narrow region of constant DM density defined by $${M_{\rm DM}}\propto R^3_{\rm 0, DM}$$. Moreover, the likelihood distribution of the initial stellar size R0, * and MBH has two major features, which correspond to two major class of allowed configurations. First, there is a vertical region with log (R0, *) close to the average size of the observed dwarfs with log (MBH) close to 0. This is a class of models with little heating and, hence, insignificant size evolution for the stars. Secondly, there is a horizontal band where log (MBH) is close to 1 and log (R0, *) close to 1. This is the allowed parameter space where a substantial heating by PBHs expands the initially more compact stellar distribution up to its present size. From these distributions, the most probable results of the four parameters are: $$\rm {log R_{\rm 0, DM}} = 2.77\, \rm {M_{\rm {\odot }}}$$, logR0, * = 1.47 pc, $$\rm {log {M_{\rm DM}}} = 9.12\, \rm {M_{\rm {\odot }}}$$, and $$\rm {log {M_{\rm BH}}} = 0.93\, \rm {M_{\rm {\odot }}}$$. Similarly, constraints using a Dehnen halo with a cusp (γ = 1) at t = 0, as shown in Fig. 6, are: $$\rm {log R_{\rm 0, DM}} = 3.34\, \rm {M_{\rm {\odot }}}$$, logR0, * = 1.46 pc, $$\rm {log {M_{\rm DM}}} = 8.92\, \rm {M_{\rm {\odot }}}$$, and $$\rm {log {M_{\rm BH}}} = 0.78\, \rm {M_{\rm {\odot }}}$$, and those with an NFW profile, as shown in Fig. 7, are: $$\rm {log R_{\rm 0, DM}} = 3.15\, \rm {M_{\rm {\odot }}}$$, logR0, * = 1.46 pc, $$\rm {log {M_{\rm DM}}} = 9.23\, \rm {M_{\rm {\odot }}}$$, and $$\rm {log {M_{\rm BH}}} = 0.76\, \rm {M_{\rm {\odot }}}$$. These results are listed in Table 2. Figure 6. View largeDownload slide Same as Fig. 5, but for a Dehnen halo with a cusp (γ = 1) at t = 0. The prior distribution of log (R0, DM) is rescaled to $${\sim } \mathcal {U}(\log (200), \log (10\,000))$$ in order to cover the allowed PBH halo parameter space. Figure 6. View largeDownload slide Same as Fig. 5, but for a Dehnen halo with a cusp (γ = 1) at t = 0. The prior distribution of log (R0, DM) is rescaled to $${\sim } \mathcal {U}(\log (200), \log (10\,000))$$ in order to cover the allowed PBH halo parameter space. Figure 7. View largeDownload slide Same as Fig. 6, but for a NFW halo truncated at 10 kpc at t = 0. Figure 7. View largeDownload slide Same as Fig. 6, but for a NFW halo truncated at 10 kpc at t = 0. Table 2. Results from 105 MCMC realizations assuming cored and cuspy DM density profiles. Model  log(R0, DM)  log(R0, *)  log(MDM)  log(MBH)  Denen core  $$2.77^{+0.17}_{-0.22}$$  $$1.47^{+0.05}_{-0.23}$$  $$9.12^{+0.42}_{-0.61}$$  $$0.93^{+0.23}_{-0.54}$$  Denen cusp  $$3.34^{+0.26}_{-0.28}$$  $$1.46^{+0.05}_{-0.27}$$  $$8.92^{+0.55}_{-0.44}$$  $$0.78^{+0.11}_{-0.41}$$  NFW cusp  $$3.15^{+0.45}_{-0.38}$$  $$1.46^{+0.043}_{-0.278}$$  $$9.23^{+0.33}_{-0.62}$$  $$0.76^{+0.14}_{-0.38}$$  Model  log(R0, DM)  log(R0, *)  log(MDM)  log(MBH)  Denen core  $$2.77^{+0.17}_{-0.22}$$  $$1.47^{+0.05}_{-0.23}$$  $$9.12^{+0.42}_{-0.61}$$  $$0.93^{+0.23}_{-0.54}$$  Denen cusp  $$3.34^{+0.26}_{-0.28}$$  $$1.46^{+0.05}_{-0.27}$$  $$8.92^{+0.55}_{-0.44}$$  $$0.78^{+0.11}_{-0.41}$$  NFW cusp  $$3.15^{+0.45}_{-0.38}$$  $$1.46^{+0.043}_{-0.278}$$  $$9.23^{+0.33}_{-0.62}$$  $$0.76^{+0.14}_{-0.38}$$  View Large A comparison of the constraints using the three different halo profiles shows that the results are very close to each other, and that these galaxies are DM-dominated: the posterior distribution of the central DM density lies in a narrow range of 1–2 M⊙pc − 3 while the stellar density is two orders of magnitude lower. We also experimented with models having no DM component at all. In this case, the velocity dispersion of stars would be much lower than the observed value for the given combination of total stellar mass and half-mass radius. By contrast, globular clusters have a far larger stellar mass while occupying the same region in the R*–σ plane. The existence of DM also ensures the logic of using ultra-faint dwarf galaxies to constrain PBH mass is self-consistent. The constraint on PBH mass in the massive end comes from a maximum heating rate, which is more stringent than that from the lower mass end, mainly due to the assumptions upon which the stellar component was initialized, as we discuss in Section 4. This can also be seen in the contour of R*, 0 and log (MBH), which is a continuous distribution. Models with R*, 0 between 20 and 30 pc are certainly allowed, with a modest heating from PBH that brings the stellar size and velocity dispersion close to the observed values. Note that the peak of MBH does not correspond to the peak of R*, 0. Therefore, the likelihood surface in the plane of R0, * and MBH is more informative, and the constraints on MBH is limited by our prior knowledge of R0, *. The estimate can be improved if we have better prior information about the initial size of the stellar component. There is a clear dichotomy in the spread in metallicity between dwarf galaxies and globular clusters of the same luminosities (see Willman & Strader 2012, most of the observed globular clusters show <0.1 dex metallicity spread). The large metallicity spread (>0.2 dex) among the member stars in the dwarf galaxies could impose more stringent constraints on the initial stellar component size. Gas confined in a small volume would lead to efficient mixing due to multiple supernova events, hence a more uniform metallicity distribution is incompatible with observations. If the initial size R0, * is above 20 pc, then MBH would fall below 10 M⊙. 4 DISCUSSIONS 4.1 Model limitations In our models, the strongest prior is the total DM halo mass, which fulfills the requirement that these haloes should be just slightly more massive than the minimum mass required for star formation. For the other three free parameters, we have adopted rather agnostic prior distributions. The initial distribution of stellar component, which follows a Plummer sphere, is a reasonable assumption. Our model assumes that DM halo has a density core. This is a reasonable assumption given that the collisional effect on a cusped density profile is well understood. To demonstrate this, we have run two additional sets of MCMC realizations with a Dehnen density cusp (γ = 1) at t = 0 and an NFW cuspy profile. As listed in Table 2 and illustrated in Figs 5–7, the final distributions of the four parameters are very close among all three models, although the peak of log (MBH) drops slightly from 0.93 to 0.76. This is due to extra heating as the central velocity dispersion is low for a cuspy density profile. Our choice of prior distribution of total DM mass is also a reasonable one. As illustrated in Fig. 5, the total DM halo mass MDM does not have strong correlation with either MBH or R*, 0. The covariance between MDM and R0, DM covers a region defined by a constant DM density in the range of 1–2 M⊙pc − 3. In contrast, PBH mass is strongly correlated with the initial size of the stellar component. Therefore, our constraints on PBH mass is unaffected by the prior distribution of total DM mass. Moreover, if more compact ultra-faint dwarfs are to be discovered in the near future, their cumulative number could impose yet another constraint against the possibility of massive M⊙ PBHs as DM. There are several limitations in our approach that need to be stressed. First, we treat the haloes in isolation while in reality they are in highly tumultuous galactic environments, subjected to various mass loss processes (Zhu et al. 2016). Secondly, we use the measurements of the compact ultra-faint dwarfs as an input for the likelihood function and treat the reported uncertainties as the widths of Gaussians in the likelihood function. We have used a smaller half-mass size for Ret II and Hor I between those reported by Koposov et al. (2015) and Bechtol et al. (2015). Accurate determination of the membership of stars in such faint galaxies is an extremely challenging task. In addition, unresolved binary stars may significantly inflate the observed velocity dispersion in these cold systems, although they are unlikely to completely dominate the measurements (McConnachie & Côté 2010). Should the true size of these stellar systems change, our modelling and conclusions need to be revisited accordingly. In our model, there is a built-in bias which puts more emphasis on the more massive PBH for the following reason. If the collisional relaxation effect is negligible, only a limited region around rh = 40 pc and σ = 4 km s−1 is allowed. On the other hand, if the model starts with a smaller initial size (with a smaller velocity dispersion), a modest heating will bring both the size and the velocity dispersion closer to the observed values. Therefore, some heating effect can be slightly more favoured than no heating effect if there are more models with modest heating allowed than those with little heating. 4.2 Comparison with previous works Our analysis improves upon the previous calculations by Brandt (2016) and Koushiappas & Loeb (2017) by using a FP code and an MCMC approach. In particular, the FP code self-consistently evolved the two components system with two mass species. In the previous studies, the velocity dispersion profiles of the two components are assumed to be the same, while our modelling with the FP code clearly shows this is not the case. Moreover, our calculations show that the heating rate from the massive species on the light species is slightly less efficient than t0.5, a result of partial equipartition of energy due to the dominant contribution from PBHs. In addition, the combination of the FP code and an MCMC approach is able to utilize both size and kinematics of the observed ultra-compact dwarf galaxies by efficiently sampling the parameter space. We note that the stellar density profile by the end of the integration contains a ‘diffused’ core in our FP calculation. This is not the same as the ring profile predicted by Koushiappas & Loeb (2017). In Fig. 8, we demonstrate the effect of PBH mass on the evolution of stellar surface density. The stars initially distribute in a Plummer sphere with a mass of M* = 350 M⊙ and a scale radius of R0, * = 20 pc. The evolved projected surface density profiles of stars at t = 12 Gyr are shown for different PBH mass ranges from MBH = 3–100 M⊙. As the PBH mass increases, the two-body heating intensifies, which leads to a stronger evolution in the core, featuring a larger core size and a lower surface density in the central region while an increase in the outer part. Clearly, the size evolution in our simulations is global, while Koushiappas & Loeb (2017) assumes that each shell evolved independently. In addition, we use an isotropic velocity distribution, which is different from the highly circular orbits, an implicit assumption made in Koushiappas & Loeb (2017). These different assumptions and treatments result in different stellar surface density profiles between our simulations and the previous analytical calculations. Figure 8. View largeDownload slide Projected stellar surface density profile from the simulation at t = 12 Gyr with different PBH masses, as represented by different colours. The initial stellar population is set up as a Plummer sphere with M* = 350 M⊙ and R0, * = 20 pc. For $$\rm {M_{BH} = 3\, \rm {M_{\rm {\odot }}}}$$, there is little evolution of the density profile for stars. As the PBH mass increases, there is a stronger evolution in the core, resulting in a larger core size and a lower surface density in the central region. Figure 8. View largeDownload slide Projected stellar surface density profile from the simulation at t = 12 Gyr with different PBH masses, as represented by different colours. The initial stellar population is set up as a Plummer sphere with M* = 350 M⊙ and R0, * = 20 pc. For $$\rm {M_{BH} = 3\, \rm {M_{\rm {\odot }}}}$$, there is little evolution of the density profile for stars. As the PBH mass increases, there is a stronger evolution in the core, resulting in a larger core size and a lower surface density in the central region. Our model assumes that the DM is completely composed of PBHs. It was suggested by Brandt (2016) that PBH above 10 M⊙ may be firmly excluded as the primary candidate of DM in dwarf galaxies. The constraint from Koushiappas & Loeb (2017) is slightly stronger, with masses greater than 6 M⊙ excluded. The constraint derived from our analysis is slightly less restrictive compared to these two studies. 4.3 Observational signatures of PBHs with millisecond time delay Our analysis is consistent with the suggestion that PBHs in the mass range of 25–100 M⊙ are ruled out as the main component of DM (Brandt 2016; Carr et al. 2017). If the LIGO detections are from mergers of PBH binaries of ∼ 30 M⊙, it is only possible if they are in the massive tail of an extended PBH mass function (Magee & Hanna 2017). However, any PBHs with a mass substantially higher than 10 M⊙ would violate the constraints from the dynamical analysis of dwarf galaxies. In the near future, aLIGO observations may soon provide constraints on PBHs around 10 M⊙. It remains possible that some fraction of the DM may be PBHs in this mass range. Apart from the contribution to the total DM budget, the existence of PBHs will provide important constraints on certain inflation theories. In addition, gravitational lensing due to ∼10 M⊙ PBHs will induce a time delay on the order of milliseconds (Mao 1992; Muñoz et al. 2016). There are at least two known astrophysical phenomena that can produce detectable signals on this rapid time-scale. One is fast radio burst (Lorimer et al. 2007). If fast radio bursts have indeed cosmological origin (Spitler et al. 2016), then one would expect repeated bursts with millisecond delay from the same location of the sky within a large number (Fialkov & Loeb 2017) of bursts due to the PBH lenses (Muñoz et al. 2016). Another promising source in the kHz regime for the detection of PBHs would be gravitational waves from stellar mass BH mergers. Similar to the mechanism of femtolensing of γ-ray burst (Gould 1992), millisecond delays could induce a characteristic interference pattern in the detected waveforms in the merger phase due to the interference of two lensed images. Perhaps gravitational waves are better sources than radio bursts since they do not suffer from interstellar dispersions. This technique could be feasible with the next generation detectors such as the ‘Cosmic Explorer’ (Abbott et al. 2017). With a large number of gravitational wave detections, we could place a strong constraint on the mass fraction of PBHs in ΩM, since the lensing optical depth is at the same order as the mass fraction of PBHs for cosmological sources (Press & Gunn 1973). 5 SUMMARY In this study, we improve the recent analytical work of Brandt (2016) and Koushiappas & Loeb (2017) on the investigation of the possibility of DM consisting entirely of PBHs, by combining Fokker–Planck simulations with MCMC realizations. The former method, implemented in the code phaeflow, accurately follows the joint evolution of stars and PBH-DM driven by two-body relaxation, for a given choice of parameters. By combining it with the MCMC approach, we explore the parameter space and use the half-light radius and central stellar velocity dispersions of observed compact ultra-faint dwarf galaxies as joint priors in the MCMC sampling to constrain the PBH-DM parameters. Our findings can be summarized as follows: For a DM halo in the mass range of $$10^6\text{--}10^9\, \rm {M_{\rm {\odot }}}$$, if the DM is composed of PBHs entirely, collisional relaxation quickly transforms a cusp to a cored density profile. For a system with stars and PBH-DM, two-body relaxation between the stars and PBH-DM drives up the stellar core size, and heating from the PBHs increases the central stellar velocity dispersion. The size evolution of the stellar core has a growth rate of rh ∝ t0.4, which is slightly slower than that from the full energy equipartition (rh ∝ t0.5). The joint constraints from the observed half-light radius and central stellar velocity dispersions of compact ultra-faint dwarf galaxies leave a narrow parameter space of DM haloes consist of PBHs. Results from the MCMC simulations suggest that the central density of DM halo has a narrow range of 1–2 M⊙pc − 3, and the mass of PBHs is most likely within 2–14 M⊙, if the DM is completely made of PBH. We note that our constraint of PBHs is slightly less restrictive than those of Brandt (2016) and Koushiappas & Loeb (2017). Although PBH with mass above 10 M⊙ can be excluded as the primary candidate of DM from our results, we emphasize that our models are inconclusive, and that the nature of DM remains elusive. ACKNOWLEDGEMENTS We thank the anonymous referee for a constructive report which has helped improve our paper. It is our great pleasure to thank Avi Loeb, Lars Hernquist, Savvas Koushiappas and Peter Mészáros for stimulating discussions. YL acknowledges support from NSF grants AST-0965694, AST-1009867, AST-1412719, and MRI-1626251. The numerical computations and data analysis in this paper have been carried out on the CyberLAMP cluster supported by MRI-1626251, operated and maintained by the Institute for CyberScience at the Pennsylvania State University, as well as the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. The Institute for Gravitation and the Cosmos is supported by the Eberly College of Science and the Office of the Senior Vice President for Research at the Pennsylvania State University. EV acknowledges support from the European Research Council under the 7th Framework Programme (grant 321067). Footnotes 1 phaseflow is part of the AGAMA library for galaxy modelling, available from https://github.com/GalacticDynamics-Oxford/Agama/. REFERENCES Abbott B. P. et al.  , 2017, Class. Quantum Gravity , 34, 044001 https://doi.org/10.1088/1361-6382/aa51f4 CrossRef Search ADS   Afshordi N., McDonald P., Spergel D. N., 2003, ApJ , 594, L71 https://doi.org/10.1086/378763 CrossRef Search ADS   Amorisco N. C., 2017, AJ , 844, 16 CrossRef Search ADS   Bechtol K. et al.  , 2015, ApJ , 807, 50 https://doi.org/10.1088/0004-637X/807/1/50 CrossRef Search ADS   Bianchini P., van de Ven G., Norris M. A., Schinnerer E., Varri A. L., 2016, MNRAS , 458, 3644 https://doi.org/10.1093/mnras/stw552 CrossRef Search ADS   Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition . Princeton Univ. Press, Princeton, NJ Bird S., Cholis I., Muñoz J. B., Ali-Haïmoud Y., Kamionkowski M., Kovetz E. D., Raccanelli A., Riess A. G., 2016, Phys. Rev. Lett. , 116, 201301 https://doi.org/10.1103/PhysRevLett.116.201301 CrossRef Search ADS PubMed  Brandt T. D., 2016, ApJ , 824, L31 https://doi.org/10.3847/2041-8205/824/2/L31 CrossRef Search ADS   Carr B. J., Hawking S. W., 1974, MNRAS , 168, 399 https://doi.org/10.1093/mnras/168.2.399 CrossRef Search ADS   Carr B., Kühnel F., Sandstad M., 2016, Phys. Rev. D , 94, 083504 https://doi.org/10.1103/PhysRevD.94.083504 CrossRef Search ADS   Carr B., Raidal M., Tenkanen T., Vaskonen V., Veermäe H., 2017, Phys. Rev. D , 96, 023514 https://doi.org/10.1103/PhysRevD.96.023514 CrossRef Search ADS   Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS , 454, 2981 https://doi.org/10.1093/mnras/stv2165 CrossRef Search ADS   Contenta F. et al.  , 2017, MNRAS , preprint (arXiv:1705.01820) Dehnen W., 1993, MNRAS , 265, 250 https://doi.org/10.1093/mnras/265.1.250 CrossRef Search ADS   Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Stinson G. S., Knebe A., 2014, MNRAS , 441, 2986 https://doi.org/10.1093/mnras/stu729 CrossRef Search ADS   Diemer B., Kravtsov A. V., 2015, ApJ , 799, 108 https://doi.org/10.1088/0004-637X/799/1/108 CrossRef Search ADS   Fialkov A., Loeb A., 2017, ApJ , 846, L27 https://doi.org/10.3847/2041-8213/aa8905 CrossRef Search ADS   Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP , 125, 306 https://doi.org/10.1086/670067 CrossRef Search ADS   Frampton P. H., Kawasaki M., Takahashi F., Yanagida T. T., 2010, J. Cosmol. Astropart. Phys. , 4, 023 https://doi.org/10.1088/1475-7516/2010/04/023 CrossRef Search ADS   Garrison-Kimmel S., Bullock J. S., Boylan-Kolchin M., Bardwell E., 2017, MNRAS , 464, 3108 https://doi.org/10.1093/mnras/stw2564 CrossRef Search ADS   Gnedin N. Y., Kaurov A. A., 2014, ApJ , 793, 30 https://doi.org/10.1088/0004-637X/793/1/30 CrossRef Search ADS   Gould A., 1992, ApJ , 386, L5 https://doi.org/10.1086/186279 CrossRef Search ADS   Hawking S., 1971, MNRAS , 152, 75 https://doi.org/10.1093/mnras/152.1.75 CrossRef Search ADS   Hernquist L., 1990, ApJ , 356, 359 https://doi.org/10.1086/168845 CrossRef Search ADS   Hinton S., 2016, JOSS , 1, 45 CrossRef Search ADS   Hirano S., Sullivan J. M., Bromm V., 2017, MNRAS , 473, L6 CrossRef Search ADS   Hu W., Barkana R., Gruzinov A., 2000, Phys. Rev. Lett. , 85, 1158 https://doi.org/10.1103/PhysRevLett.85.1158 CrossRef Search ADS PubMed  Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, Phys. Rev. D , 95, 043541 https://doi.org/10.1103/PhysRevD.95.043541 CrossRef Search ADS   Kashlinsky A., 2016, ApJ , 823, L25 https://doi.org/10.3847/2041-8205/823/2/L25 CrossRef Search ADS   Khlopov M. Y., 2010, Res. Astron. Astrophys. , 10, 495 https://doi.org/10.1088/1674-4527/10/6/001 CrossRef Search ADS   Koposov S. E., Belokurov V., Torrealba G., Evans N. W., 2015, ApJ , 805, 130 https://doi.org/10.1088/0004-637X/805/2/130 CrossRef Search ADS   Koushiappas S. M., Loeb A., 2017, Phys. Rev. Lett. , 119, 041102 https://doi.org/10.1103/PhysRevLett.119.041102 CrossRef Search ADS PubMed  Lorimer D. R., Bailes M., McLaughlin M. A., Narkevic D. J., Crawford F., 2007, Science , 318, 777 https://doi.org/10.1126/science.1147532 CrossRef Search ADS PubMed  Ma X. et al.  , 2017, MNRAS , preprint (arXiv:1706.06605) Magee R., Hanna C., 2017, AJ , 845, 5 CrossRef Search ADS   Mao S., 1992, ApJ , 389, L41 https://doi.org/10.1086/186344 CrossRef Search ADS   McConnachie A. W., Côté P., 2010, ApJ , 722, L209 https://doi.org/10.1088/2041-8205/722/2/L209 CrossRef Search ADS   Meszaros P., 1975, A&A , 38, 5 Muñoz J. B., Kovetz E. D., Dai L., Kamionkowski M., 2016, Phys. Rev. Lett. , 117, 091301 https://doi.org/10.1103/PhysRevLett.117.091301 CrossRef Search ADS PubMed  Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ , 490, 493 https://doi.org/10.1086/304888 CrossRef Search ADS   Okamoto T., Gao L., Theuns T., 2008, MNRAS , 390, 920 https://doi.org/10.1111/j.1365-2966.2008.13830.x CrossRef Search ADS   Plummer H. C., 1911, MNRAS , 71, 460 https://doi.org/10.1093/mnras/71.5.460 CrossRef Search ADS   Pontzen A., Governato F., 2014, Nature , 506, 171 https://doi.org/10.1038/nature12953 CrossRef Search ADS PubMed  Press W. H., Gunn J. E., 1973, ApJ , 185, 397 https://doi.org/10.1086/152430 CrossRef Search ADS   Quinlan G. D., 1996, New A , 1, 255 https://doi.org/10.1016/S1384-1076(96)00018-8 CrossRef Search ADS   Read J. I., Agertz O., Collins M. L. M., 2016, MNRAS , 459, 2573 https://doi.org/10.1093/mnras/stw713 CrossRef Search ADS   Sawala T. et al.  , 2015, MNRAS , 448, 2941 https://doi.org/10.1093/mnras/stu2753 CrossRef Search ADS   Shapiro P. R., Iliev I. T., Raga A. C., 1999, MNRAS , 307, 203 https://doi.org/10.1046/j.1365-8711.1999.02609.x CrossRef Search ADS   Spitler L. G. et al.  , 2016, Nature , 531, 202 https://doi.org/10.1038/nature17168 CrossRef Search ADS PubMed  Spitzer L., Jr, 1969, ApJ , 158, L139 https://doi.org/10.1086/180451 CrossRef Search ADS   The LIGO Scientific Collaboration, 2016a, Phys. Rev. Lett. , 116, 061102 https://doi.org/10.1103/PhysRevLett.116.061102 CrossRef Search ADS   The LIGO Scientific Collaboration, 2016b, Phys. Rev. Lett. , 116, 241103 https://doi.org/10.1103/PhysRevLett.116.241103 CrossRef Search ADS   The LIGO Scientific Collaboration, 2017a, Phys. Rev. Lett. , 119, 141101 CrossRef Search ADS   The LIGO Scientific Collaboration, 2017b, Phys. Rev. Lett. , 118, 221101 https://doi.org/10.1103/PhysRevLett.118.221101 CrossRef Search ADS   Trenti M., van der Marel R., 2013, MNRAS , 435, 3272 https://doi.org/10.1093/mnras/stt1521 CrossRef Search ADS   Vasiliev E., 2017, ApJ , 848, 10 https://doi.org/10.3847/1538-4357/aa8cc8 CrossRef Search ADS   Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS , 423, 3740 https://doi.org/10.1111/j.1365-2966.2012.21182.x CrossRef Search ADS   Vogelsberger M., Zavala J., Cyr-Racine F.-Y., Pfrommer C., Bringmann T., Sigurdson K., 2016, MNRAS , 460, 1399 https://doi.org/10.1093/mnras/stw1076 CrossRef Search ADS   Weinberg D. H., Bullock J. S., Governato F., Kuzio de Naray R., Peter A. H. G., 2015, Proc. Natl. Acad. Sci. , 112, 12249 https://doi.org/10.1073/pnas.1308716112 CrossRef Search ADS   Wheeler C., Oñorbe J., Bullock J. S., Boylan-Kolchin M., Elbert O. D., Garrison-Kimmel S., Hopkins P. F., Kereš D., 2015, MNRAS , 453, 1305 https://doi.org/10.1093/mnras/stv1691 CrossRef Search ADS   Willman B., Strader J., 2012, AJ , 144, 76 https://doi.org/10.1088/0004-6256/144/3/76 CrossRef Search ADS   Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., Muñoz R. R., Simon J. D., Avedo F. F., 2010, MNRAS , 406, 1220 Yoshida N., Sokasian A., Hernquist L., Springel V., 2003, ApJ , 591, L1 https://doi.org/10.1086/376963 CrossRef Search ADS   Zhao H., 1996, MNRAS , 278, 488 https://doi.org/10.1093/mnras/278.2.488 CrossRef Search ADS   Zhu Q., Marinacci F., Maji M., Li Y., Springel V., Hernquist L., 2016, MNRAS , 458, 1559 https://doi.org/10.1093/mnras/stw374 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Primordial black holes as dark matter: constraints from compact ultra-faint dwarfs

, Volume 476 (1) – May 1, 2018
10 pages

/lp/ou_press/primordial-black-holes-as-dark-matter-constraints-from-compact-ultra-hErpISTBid
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty079
Publisher site
See Article on Publisher Site

### Abstract

Abstract The ground-breaking detections of gravitational waves from black hole mergers by LIGO have rekindled interest in primordial black holes (PBHs) and the possibility of dark matter being composed of PBHs. It has been suggested that PBHs of tens of solar masses could serve as dark matter candidates. Recent analytical studies demonstrated that compact ultra-faint dwarf galaxies can serve as a sensitive test for the PBH dark matter hypothesis, since stars in such a halo-dominated system would be heated by the more massive PBHs, their present-day distribution can provide strong constraints on PBH mass. In this study, we further explore this scenario with more detailed calculations, using a combination of dynamical simulations and Bayesian inference methods. The joint evolution of stars and PBH dark matter is followed with a Fokker–Planck code phaseflow. We run a large suite of such simulations for different dark matter parameters, then use a Markov chain Monte Carlo approach to constrain the PBH properties with observations of ultra-faint galaxies. We find that two-body relaxation between the stars and PBH drives up the stellar core size, and increases the central stellar velocity dispersion. Using the observed half-light radius and velocity dispersion of stars in the compact ultra-faint dwarf galaxies as joint constraints, we infer that these dwarfs may have a cored dark matter halo with the central density in the range of 1–2 M⊙pc − 3, and that the PBHs may have a mass range of 2–14 M⊙ if they constitute all or a substantial fraction of the dark matter. methods: numerical, dark matter 1 INTRODUCTION In the standard lambda cold dark matter (ΛCDM) cosmology, the nature of the dark matter (DM) remains elusive. The possibility of DM being composed of primordial black holes (PBHs; Hawking 1971), arising naturally from certain inflation models (e.g. Frampton et al. 2010), has been investigated in the context of galaxy formation early on (e.g. Carr & Hawking 1974; Meszaros 1975), and it has not yet been completely ruled out (see recent reviews by Khlopov 2010; Carr, Kühnel & Sandstad 2016). Recently, it was suggested that three mass windows, around 5 × 10 − 16 M⊙, 2 × 10 − 14 M⊙, and 25–100 M⊙, respectively, are still open for DM consists of PBHs (Carr et al. 2017). The recent detections of gravitational waves by LIGO are believed to originate from merging black holes with masses of order 10–50 M⊙, significantly higher than known BHs in Galactic X-ray binaries (The LIGO Scientific Collaboration et al. 2016a,b, 2017a,b). These discoveries have renewed theoretical interest in PBHs as a DM candidate. It was suggested by Bird et al. (2016) that PBHs with a mass around 30 M⊙ are loosely constrained by micro-lensing experiments and that the merger rate of binary PBHs is consistent with the LIGO detections. If true, PBHs in this mass range could be an attractive DM candidate to solve some well-known problems on the galactic scale (see Weinberg et al. 2015), in particular the ‘cusp versus core problem’ (hereafter, we refer to dark matter made of primordial black holes as PBH-DM). The PBHs would produce DM cores instead of cusps in the central galaxy because any ‘temperature inversion’ present in a cuspy profile will lead to a removal of the density cusp by increasing the velocity dispersion (Quinlan 1996). Solutions to the ‘cusp versus core problem’ based on baryonic heating of DM are shown to be mass-dependent (Di Cintio et al. 2014; Pontzen & Governato 2014; Chan et al. 2015). The survival of a globular cluster in the ultra-faint dwarf galaxy Eridanus II, however, hints to a DM core in a galaxy halo where any baryonic effect would be weak (Amorisco 2017; Contenta et al. 2017). In addition, a cosmology model with PBHs as DM indicates a rather different evolution picture from that of other types of DM, such as warm dark matter, fuzzy dark matter (Hu, Barkana & Gruzinov 2000; Hui et al. 2017), and self-interacting dark matter (Vogelsberger, Zavala & Loeb 2012). Contrary to a common feature of suppressed small-scale perturbations in these models (Vogelsberger et al. 2016), fluctuations in the number density of PBH impose additional small-scale power (Afshordi, McDonald & Spergel 2003). As a result, low-mass PBH-DM haloes will form ahead of ΛCDM haloes, while in warm/fuzzy dark matter, the DM halo formation is considerably delayed, which has a strong consequence on the onset of reionization (Yoshida et al. 2003; Hirano, Sullivan & Bromm 2017). In addition to the increased collapsed halo at high redshift, Kashlinsky (2016) has shown that gas accretion on to the PBHs could explain the observed high coherence between cosmic infrared and the soft-X-ray background. Despite these appealing features, it remains unclear what fraction of the DM can PBHs constitute, and what mass range can these PBHs have. Fortunately, a population of compact (with projected half-light radius of ∼30 pc), ultra-faint (∼1000 L⊙) dwarf galaxies has recently been discovered (e.g. Koposov et al. 2015; Bechtol et al. 2015), and they are expected to provide stringent constraints on the PBH-DM. Using analytical approaches, Brandt (2016) suggested that two-body relaxation between stars and PBH-DM would drive a substantial size evolution in galaxies, which is inconsistent with the observed small size of the ultra-faint dwarfs, while Koushiappas & Loeb (2017) argued that mass segregation would rearrange the stars into an otherwise detectable ring in the surface brightness profile of Segue I. In order to improve these analytical studies, we combine two methods, Fokker–Planck (FP) simulations and Markov chain Monte Carlo (MCMC) modelling, to study the evolution of galaxy haloes containing both stars and PBH-DM and to compare with observations. We first use an FP code to perform a large set of galaxy simulations with different parameters for the PBH-DM to follow the interaction and evolution of the stars and PBHs, then we employ an MCMC method to constrain the PBH properties using the observed half-light radius and velocity dispersion of stars in compact ultra-faint galaxies. Our modelling presents a significant improvement over the previous studies based on an analytical form of energy diffusion coefficient and other simplifying assumptions (Brandt 2016; Koushiappas & Loeb 2017). This paper is organized as follows. In Section 2, we describe the numerical code and MCMC algorithm used in our simulations and analysis. In Section 3, we present the dynamical evolution of the stars and PBHs in the halo systems, and constraints from observations. We discuss the implications of our study and its limitations in Section 4, and summarize our main findings in Section 5. 2 METHODS 2.1 Numerical code for dynamical simulations In this study, we use an FP code, phaseflow1 (Vasiliev 2017) to follow the dynamical evolution of stars and PBH-DM in the galaxy haloes. phaseflow solves the one-dimensional FP equation for the distribution function in energy, using a high-accuracy finite element method. The code is well tested with problems such as core collapse and the formation of Bahcall–Wolf cusp around a central massive black hole. phaseflow can handle multiple mass components, which enables our calculations of the evolution of a halo of stars and DM consisting of PBHs. Once the initial conditions of the two components are specified, phaseflow evolves the distribution function of all components while at the same time updating the gravitational potential. System diagnosis such as density profile and velocity dispersions are computed automatically, which greatly simplifies our analysis. Throughout the paper, we set the Coulomb logarithm (Binney & Tremaine 2008) ln Λ to be 15. The stars in the DM halo are modelled with a Plummer sphere (Plummer 1911):   \begin{eqnarray} \rho (r) = \frac{3M_{*}}{4\pi R_{0,*}^3} \bigg (1+\frac{r^2}{R_{0,*}^2} \bigg )^{-5/2}, \end{eqnarray} (1)where the stellar system is characterized by the total stellar mass M* and its scale radius R0, *, respectively. In this calculation, we set M* = 103 M⊙, and each star particle has an equal mass of 1 M⊙. The extended PBH-DM halo follows a Dehnen sphere (Dehnen 1993; Contenta et al. 2017):   \begin{eqnarray} \rho (r) = \frac{(3-\gamma ){M_{\rm DM}}}{4\pi R_{\rm 0, DM}^3} \bigg (\frac{r}{R_{\rm 0, DM}}\bigg )^{-\gamma } \bigg (1 + \frac{r}{R_{\rm 0, DM}}\bigg )^{\gamma -4}, \end{eqnarray} (2)where MDM is the total mass and R0, DM the scale radius of the halo, respectively. The Dehnen sphere is a special case of the more general α–β–γ density profile (Zhao 1996), where the density distribution is determined by the inner slope γ and the outer slope β and the steepness parameter α, which control how quickly the inner slope transits into the outer slope. Note that when γ = 1, the density distribution becomes the Hernquist profile (Hernquist 1990). In order to compare with previous studies, we also consider the Navarro–Frenk–White (NFW) profile (Navarro, Frenk & White 1997) in our simulations. 2.2 Markov chain Monte Carlo for parameter constraints 2.2.1 Observational data and the likelihood function The MCMC is an effective method to explore multiple key parameters in a complex system. In our model, we use both the stellar half-light radius and the stellar velocity dispersion of the compact ultra-faint dwarfs to constrain the parameter space, which requires ∼105 MCMC realizations. Fortunately, phaseflow is an extremely efficient code, and it takes only several seconds to run a single model. Thus, we combine a large grid of FP simulations performed using phaseflow and MCMC calculations using the emcee code. The emcee code uses an affine-invariant ensemble sampler to efficiently random walk the parameter space (Foreman-Mackey et al. 2013). The likelihood function of our modelling is described in the following form:   \begin{eqnarray} \mathcal {L} = \prod _{i=1}^{N} \mathcal {N} (r_{\rm h, M}, r_{\rm h, O}, \epsilon _{\rm h, O})\, \mathcal {N} (\sigma _{\rm M}, \sigma _{\rm O}, \epsilon _{\rm \sigma , O}), \end{eqnarray} (3)where $$\mathcal {N}$$ is a normal distribution, N is the sample size, rh, M is the model output of 3D half-mass radius, rh, O, and σO are the observed half-light radius and velocity dispersion, respectively. The associated uncertainties are εh, O and εσ, O. We use the top five ultra-faint dwarf galaxies compiled by Brandt (2016), which have measurements of both size and stellar velocity dispersion available. The data are repeated here in Table 1 for convenience but we refer the readers to Brandt (2016) for details and references of each measurement. Note that there are discrepancies in the reported half-mass radius of the two newly discovered Ret II and Hor I between Koposov et al. (2015) and Bechtol et al. (2015). In each case, we have chosen the smaller value of the two measurements. The reason for such a sample is based on the implicit assumption (inherited from Brandt 2016) that these five galaxies form a distinct class of faint dwarfs with similar size and velocity dispersion, which can be compared to a single theoretical model. Lastly, we apply a factor of 1.3 to de-project the observed half-mass radius to 3D following Wolf et al. (2010). Table 1. A sample of observed compact ultra-faint dwarf galaxies used in our modelling from Brandt (2016). Galaxy name  Projected rh (pc)  σ* (km s−1)  LV (L⊙)  Wil I  25 ± 6  $$4.3^{+2.3}_{-1.3}$$  1000  Seg I  $$29^{+8}_{-5}$$  $$3.9^{+0.8}_{-0.8}$$  300  Seg II  35 ± 3  $$3.4^{+2.5}_{-1.2}$$  900  Ret II  $$32^{+2}_{-1}$$  $$3.2^{+1.6}_{-0.5}$$  1500  Hor I  $$25^{+9}_{-4}$$  $$4.9^{+2.8}_{-0.9}$$  2000  Galaxy name  Projected rh (pc)  σ* (km s−1)  LV (L⊙)  Wil I  25 ± 6  $$4.3^{+2.3}_{-1.3}$$  1000  Seg I  $$29^{+8}_{-5}$$  $$3.9^{+0.8}_{-0.8}$$  300  Seg II  35 ± 3  $$3.4^{+2.5}_{-1.2}$$  900  Ret II  $$32^{+2}_{-1}$$  $$3.2^{+1.6}_{-0.5}$$  1500  Hor I  $$25^{+9}_{-4}$$  $$4.9^{+2.8}_{-0.9}$$  2000  View Large 2.2.2 Selections of prior distributions Total dark matter mass MDM: Compact ultra-faint dwarf galaxies from our sample have the least amount of stars (∼1000 L⊙) ever known, but their total mass remains unknown. It was suggested that they may reside in the least massive haloes where atomic cooling is efficient, and that their total mass is around 109 M⊙ (Sawala et al. 2015; Wheeler et al. 2015; Zhu et al. 2016), since such a halo mass is consistent with the M*–MDM relations at the low mass end (Garrison-Kimmel et al. 2017). Recently, Ma et al. (2017) resolved the formation of such galaxies in their simulations, and reported that the least massive halo, which has a mass of 108M⊙ at z = 5, had already stopped star formation since z = 8. On the other hand, Read, Agertz & Collins (2016) modelled isolated dwarf galaxies and found that haloes with a total mass of 108M⊙ were able to match several dwarf satellites of the Milky Way. Therefore, we choose a weakly informative prior of MDM, which follows a log-normal distribution:   \begin{eqnarray} \log ({M_{\rm DM}}) \sim \mathcal {N}(9, 0.5). \end{eqnarray} (4) The choice of this prior is well justified in that the physics involved is relatively straightforward. These low-mass galaxies must be just massive enough, above the thresholds set by atomic cooling and UV background radiation, for stars to form (e.g. Okamoto, Gao & Theuns 2008; Gnedin & Kaurov 2014). The width of the distribution also ensures that haloes down to 2–3 × 108 M⊙ are covered. DM scale radius R0, DM:From the discussion below (Section 3.1), we argue that the inner density profile of the halo will form a core (γ = 0) due to relaxation effect. However, the size of the scale radius needs to be treated as another free parameter in our model. This is largely due to a lack of self-consistent collisional N-body simulations from cosmological initial conditions. The assumption that the characteristic scales of the PBH-DM haloes are similar to the ΛCDM models could lead to serious systematics. Thus, we adopt a log-uniform distribution of R*, DM for the range from 100 to 2000 pc:   \begin{eqnarray} \log (R_{\rm *, DM}) \sim \mathcal {U}(\log (100), \log (2000)). \end{eqnarray} (5) Stellar scale radius R0, *:The heating of stars due to two-body relaxation increases the scale radius of stellar distribution roughly as R* ∼ t0.5 (Brandt 2016), which indicates that the present-day size does not depend sensitively on the initial value of R0, *. We thus use a log-uniform distribution for R0, * for the range from 10 to 50 pc.   \begin{eqnarray} \log (R_{\rm 0, *}) \sim \mathcal {U}(\log (10), \log (50)). \end{eqnarray} (6) The lower limit corresponds to the largest Galactic globular clusters and the upper limit corresponds to the average size of the currently observed compact ultra-faint dwarfs. Primordial BH mass MBH: Lastly, the mass of the PBHs is not well constrained (Carr et al. 2017). PBHs above 100 M⊙ can be ruled out by CMB observation and the Galactic wide binaries (see Carr et al. 2016). However, since we are most interested in the range of [1, 100] M⊙, we assume a log-uniformly distribution for MBH in this range. 3 RESULTS 3.1 Cusp-to-core transition in PBH-DM haloes from relaxation effect In a CDM halo, the density profile usually shows a central density cusp with a slope in the range of 0 < γ < 2, and the DM velocity dispersion peaks at the scale radius. However, if the DM is composed of PBHs, two-body relaxation will soften the central cusp, and as a result of the collisional heating, the velocity dispersion in the central region will increase while the density will drop at the same time. This process leads to the rapid formation of a core (Quinlan 1996) as a result of ‘temperature inversion’, in which the colder dense cusp will be heated by a hotter envelope. Fig. 1 illustrates how this relaxation effect transforms a PBH-DM halo from a γ = 1 cusp to a γ = 0 core. The DM halo has a mass of 2 × 109 M⊙ composed of 30 M⊙ PBHs, and an initial cuspy density profile of γ = 1 and R0, DM = 500 pc. Collisional relaxation quickly removes the central density cusp and increases the central velocity dispersion. After 12 Gyr, this PBH-DM halo can be approximated by a less massive (109M⊙) one with a core density profile of γ = 0 and a characteristic scale radius of R0, DM = 160 pc. In addition, the evolution of an NFW halo with a scale radius of 300 pc and truncated at 10 kpc is shown in Fig. 1 for comparison. The scale radius is chosen to match the initial cuspy Dehnen profile. Within the central 200 pc, both the density and the velocity dispersion profiles of the NFW halo match perfectly with those of the cuspy Dehnen model. While the NFW profile leads different asymptotic behaviour towards large r, but the dynamical region we focus in this study is within the central ∼ 100 pc, where both NFW and the cuspy Dehnen profiles produce the same results. Figure 1. View largeDownload slide The transition from cusp to core (top panel) and the increase of central velocity dispersion (bottom panel) due to collisional relaxation effect in a PBH-DM halo. The DM halo has a mass of 2 × 109M⊙ consisting of 30 M⊙ PBHs. It has an initial density profile of γ = 1 (Hernquist profile) and R0, DM = 500 pc, and a diminishing velocity dispersion as r → 0, as represented by the thick black curves in both panels. The thin lines show the evolution of the density profile and the velocity dispersion from t = 0 to 12 Gyr with an interval of 1.5 Gyr. The red dashed lines shows the best matching PBH-DM halo of 109M⊙ with a core density profile of γ = 0 and R0, DM = 160 pc, and an increased central velocity dispersion. For comparison, we also show the evolution of an NFW halo, which is truncated at 10 kpc, at t = 0 and 12 Gyr, respectively (green dashed curves). Figure 1. View largeDownload slide The transition from cusp to core (top panel) and the increase of central velocity dispersion (bottom panel) due to collisional relaxation effect in a PBH-DM halo. The DM halo has a mass of 2 × 109M⊙ consisting of 30 M⊙ PBHs. It has an initial density profile of γ = 1 (Hernquist profile) and R0, DM = 500 pc, and a diminishing velocity dispersion as r → 0, as represented by the thick black curves in both panels. The thin lines show the evolution of the density profile and the velocity dispersion from t = 0 to 12 Gyr with an interval of 1.5 Gyr. The red dashed lines shows the best matching PBH-DM halo of 109M⊙ with a core density profile of γ = 0 and R0, DM = 160 pc, and an increased central velocity dispersion. For comparison, we also show the evolution of an NFW halo, which is truncated at 10 kpc, at t = 0 and 12 Gyr, respectively (green dashed curves). In order to determine the time-scale of the transition from cusp-to-core profiles, we set up cuspy DM haloes consisting of 30 M⊙ PBHs with a total mass in the range from 105 to 109 M⊙ at redshift z = 10, based on the mass-concentration relation from Diemer & Kravtsov (2015). We evolve the PBH-DM haloes using the FP code phaseflow. These simulations show that the collisional relaxation effect removes the central cusp almost instantaneously, and the DM core grows quickly in size. For instance, it takes only ∼0.05 Gyr for all the haloes to develop a core of ∼10 pc, and ∼0.18 Gyr to reach ∼20 pc. If formed via hierarchical assembly, it would be difficult for DM haloes of $$10^8\text{--}10^9\, \rm {M_{\rm {\odot }}}$$ to retain density cusps in the first place since the central density of their less massive progenitors has already been lowered. We speculate that the density profile of a PBH-DM halo, similar to self-interacting DM (Vogelsberger et al. 2012), could be described by a truncated singular sphere with a sizable core (Shapiro, Iliev & Raga 1999). In the absence of a self-consistent treatment of density profile of PBH-DM haloes, we thus use a cored density profile with γ = 0 for any initial DM scale radius R0, DM as our default choice for simulations in the following sections. In the Discussion Section, we will include a test using γ = 1 to examine the outcome with cuspy density profiles. 3.2 Stars in PBH-DM haloes In a system with two components of different masses, two-body encounters and energy exchange would lead to change in the distribution of density and energy. To investigate the evolution of a DM halo with stars and PBHs, we set up a DM halo with a mass of 2 × 109 M⊙, which contains 103 M⊙ of stars, each being 1 M⊙, and the DM is composed of 30 M⊙ PBHs. The Plummer sphere for the stars has a total mass of 103M⊙ and a scale radius of 5 pc, while the Dehnen sphere for the DM has a total mass of 2 × 109M⊙ and a scale radius R0, DM = 500 pc. The system is then followed with the FP code phaseflow, and the results are shown in Fig. 2. Due to the huge difference in both mass and size, the two components have different initialization in both density profile and velocity dispersion, and strong two-body relaxation and energy exchange leads to significant change in the stellar density profile and velocity dispersion. Because the PBH is much more massive than the star, heating from PBHs on the stars is the dominant driver for the evolution of the stellar component. After 12 Gyr, the central stellar density drops by nearly 3 orders of magnitude, as shown in the upper panel of Fig. 2, while the central stellar velocity dispersion increases by a factor of ∼8, as in the lower panel. Figure 2. View largeDownload slide Effects of two-body relaxation and energy exchange on the density profile (top panel) and velocity dispersion (bottom panel) of a two-component halo system with stars and PBH-DM. The halo has a total DM mass of 2 × 109 M⊙ composed of 30 M⊙ PBHs, and a total stellar mass of 103 M⊙ consisting of 1 M⊙ stars. The blue curves represent the PBH-DM component at t = 0 (dashed lines) and t = 12 Gyr (solid lines), while the red curves represent the stellar components at t = 0 (dashed lines) and t = 12 Gyr (solid lines). Figure 2. View largeDownload slide Effects of two-body relaxation and energy exchange on the density profile (top panel) and velocity dispersion (bottom panel) of a two-component halo system with stars and PBH-DM. The halo has a total DM mass of 2 × 109 M⊙ composed of 30 M⊙ PBHs, and a total stellar mass of 103 M⊙ consisting of 1 M⊙ stars. The blue curves represent the PBH-DM component at t = 0 (dashed lines) and t = 12 Gyr (solid lines), while the red curves represent the stellar components at t = 0 (dashed lines) and t = 12 Gyr (solid lines). As the density profile of the stars slowly diffuses out, the half-light radius increases. This process depends on the density and velocity dispersion profile of the PBH halo. We note that in Fig. 2, the stellar velocity dispersion is well below that of PBH-DM at t = 12 Gyr. Hence, equipartition of energy between the two mass species, i.e. σ(m)2 ∝ m−1, has not been achieved. As first pointed out by Spitzer (1969), full energy equipartition is possible only if the mass fraction of heavy species is below some critical value. In our setup, the total mass ratio between PBHs and stars is well above the critical value ∼0.001(= 0.16(1/30)1.5). As a result, the central regime will be devoid of stars, leading to only a ‘partial equipartition’. Recently, Trenti & van der Marel (2013) and Bianchini et al. (2016) reported partial equipartition in their direct N-body simulations of globular clusters. Of course, our example is quite extreme since the mass density profile is completely dominated by PBHs at almost all radii. The size increase due to heating of PBH appears to be slightly slower than the analytical result obtained by Brandt (2016). To further investigate this, we performed a set of FP simulations with phaseflow for different initial stellar density distributions, by varying R0, * from 3 to 15 pc, within a fixed Dehnen sphere with a total DM mass of 2 × 109M⊙ and R0, DM = 500 pc. The resulting size evolution of the different stellar components is shown in Fig. 3, in comparison with that of Brandt (2016). Our simulations show a slower growth rate, rh ∝ t0.4, than that of rh ∝ t0.5 by Brandt (2016). The heating rate of a less concentrated stellar component is also slower than those of more concentrated ones. For example, the magenta line on top of Fig. 3 shows the size evolution of a stellar core from R0, * = 15 pc, which only approaches the asymptotic t0.4 trend by the end of the integration. Our results suggest that the final size of the stellar core only depends weakly on the initial size, and that after ∼10  Gyr, a two-component PBH-DM halo of 2 × 109M⊙ would produce a stellar core of ∼50 pc, regardless of its initial size. Figure 3. View largeDownload slide The evolution of 3D half-light radius of stellar components with various initial conditions in a PBH-DM halo from our simulations, in comparison with analytical result from Brandt (2016). The different colour represents different initial condition with R0, * varying from 3, 6, 9, 12, and 15 pc, respectively. Our results show an rh ∼ t0.4 size growth rate, slower than rh ∼ t0.5 reported by Brandt (2016). Figure 3. View largeDownload slide The evolution of 3D half-light radius of stellar components with various initial conditions in a PBH-DM halo from our simulations, in comparison with analytical result from Brandt (2016). The different colour represents different initial condition with R0, * varying from 3, 6, 9, 12, and 15 pc, respectively. Our results show an rh ∼ t0.4 size growth rate, slower than rh ∼ t0.5 reported by Brandt (2016). 3.3 Parameter space of PBH-DM haloes from FP simulations In order to explore the parameter space of DM haloes of the observed compact ultra-faint dwarfs in Table 1, we use the FP code phaseflow to perform a grid of dynamical simulations of the two-component systems consisting of stars and PBH-DM with different halo mass and size. We vary the total DM mass from 108 to 5 × 109 M⊙, and vary the scale radius of the DM Dehnen sphere from 100 to 1000 pc. We consider 50 values for each of the two DM halo parameters, MDM and R0, DM, resulting in a total of 2500 simulations. We assume the DM is composed of PBHs of 30 M⊙, a total stellar mass of 103 M⊙, and an initial stellar scale radius of R0, * = 20 pc. Each system is integrated for a duration of 12 Gyr. The final half-mass radius and the projected central velocity dispersion of the stars are shown in Fig. 4, as functions of the DM parameter space of total DM masses MDM and scale radius R0, DM. The shaded regions in both the top and bottom panels indicate the allowed range of DM parameters from the observations of the compact ultra-faint dwarfs. From the stellar half-mass radius distribution in the top panel, given the same total DM mass, the spreading of the stellar component is more efficient when the scale radius is larger. While with the same scale radius, less massive haloes show more efficient heating. When the halo is less dense (either due to larger scale radius or lower mass), the initial velocity dispersions of both stars and PBHs are also much lower. In this case, the relative increase of stellar velocity dispersion is larger than in a more dense halo. The preferred range of rh ∼ 25–40 pc from observations demand a narrow range of R0, DM < 400 pc and MDM < 5 × 109 M⊙. On the other hand, from the stellar velocity dispersion in the bottom panel, the preferred range of σ ∼ 3–5 km s−1 dictates a range of 300 < R0, DM < 1000 pc and MDM > 109 M⊙. Figure 4. View largeDownload slide The distribution of 3D half-mass radius (top panel) and 1D velocity dispersion (bottom panel) of the stellar components in the PBH-DM halo parameter space, from 2500 FP simulations using the phaseflow code. There are 50 grids for each of the DM halo parameters, total DM mass MDM and scale radius R0, DM. The DM is assumed to consist of PBHs of 30 M⊙. The allowed regions by observations are highlighted and hatched in both panels which, however, have little overlap with each other. Figure 4. View largeDownload slide The distribution of 3D half-mass radius (top panel) and 1D velocity dispersion (bottom panel) of the stellar components in the PBH-DM halo parameter space, from 2500 FP simulations using the phaseflow code. There are 50 grids for each of the DM halo parameters, total DM mass MDM and scale radius R0, DM. The DM is assumed to consist of PBHs of 30 M⊙. The allowed regions by observations are highlighted and hatched in both panels which, however, have little overlap with each other. Clearly, it is difficult to reconcile the two constraints simultaneously. A smaller stellar size would require a DM halo with smaller scale radius, thus reducing the heating from relaxation by increasing the velocity dispersion. Given PBH mass of 30 M⊙ and a typical size of ∼40 pc (de-projected), there is little room in the plane of MDM and r0, DM to meet both the constraints. 3.4 MCMC constraints on PBH-DM from compact ultra-faint dwarfs In order to constrain the PBH-DM parameters from observations, we use the stellar half-mass radius and stellar velocity dispersion of the observed compact ultra-faint dwarfs in Table 1 as priors and perform a total of 105 MCMC realizations using the emcee code. The resulting posterior probability distributions of the four model parameters, the DM scale radius R0, DM, initial stellar component scale radius R0, *, total DM mass MDM, and PBH mass MBH, are shown in Fig. 5, respectively. Figure 5. View largeDownload slide Posterior probability distributions of the four PBH-DM parameters constrained by the stellar half-mass radius and stellar velocity dispersion of observed compact ultra-faint dwarfs: the DM scale radius R0, DM, initial stellar component scale radius R0, *, total DM mass MDM, and PBH mass MBH. The contours indicate the 1σ and 2σ uncertainty region. This plot is generated using the python package ChainConsumer (Hinton 2016). Figure 5. View largeDownload slide Posterior probability distributions of the four PBH-DM parameters constrained by the stellar half-mass radius and stellar velocity dispersion of observed compact ultra-faint dwarfs: the DM scale radius R0, DM, initial stellar component scale radius R0, *, total DM mass MDM, and PBH mass MBH. The contours indicate the 1σ and 2σ uncertainty region. This plot is generated using the python package ChainConsumer (Hinton 2016). The posterior surface in the plane of R0, DM and M0, DM can be understood with the help of Fig. 4, as a balance between the requirements of small sizes and small velocity dispersions. The contours in the plane of R0, DM and M0, DM roughly cover a narrow region of constant DM density defined by $${M_{\rm DM}}\propto R^3_{\rm 0, DM}$$. Moreover, the likelihood distribution of the initial stellar size R0, * and MBH has two major features, which correspond to two major class of allowed configurations. First, there is a vertical region with log (R0, *) close to the average size of the observed dwarfs with log (MBH) close to 0. This is a class of models with little heating and, hence, insignificant size evolution for the stars. Secondly, there is a horizontal band where log (MBH) is close to 1 and log (R0, *) close to 1. This is the allowed parameter space where a substantial heating by PBHs expands the initially more compact stellar distribution up to its present size. From these distributions, the most probable results of the four parameters are: $$\rm {log R_{\rm 0, DM}} = 2.77\, \rm {M_{\rm {\odot }}}$$, logR0, * = 1.47 pc, $$\rm {log {M_{\rm DM}}} = 9.12\, \rm {M_{\rm {\odot }}}$$, and $$\rm {log {M_{\rm BH}}} = 0.93\, \rm {M_{\rm {\odot }}}$$. Similarly, constraints using a Dehnen halo with a cusp (γ = 1) at t = 0, as shown in Fig. 6, are: $$\rm {log R_{\rm 0, DM}} = 3.34\, \rm {M_{\rm {\odot }}}$$, logR0, * = 1.46 pc, $$\rm {log {M_{\rm DM}}} = 8.92\, \rm {M_{\rm {\odot }}}$$, and $$\rm {log {M_{\rm BH}}} = 0.78\, \rm {M_{\rm {\odot }}}$$, and those with an NFW profile, as shown in Fig. 7, are: $$\rm {log R_{\rm 0, DM}} = 3.15\, \rm {M_{\rm {\odot }}}$$, logR0, * = 1.46 pc, $$\rm {log {M_{\rm DM}}} = 9.23\, \rm {M_{\rm {\odot }}}$$, and $$\rm {log {M_{\rm BH}}} = 0.76\, \rm {M_{\rm {\odot }}}$$. These results are listed in Table 2. Figure 6. View largeDownload slide Same as Fig. 5, but for a Dehnen halo with a cusp (γ = 1) at t = 0. The prior distribution of log (R0, DM) is rescaled to $${\sim } \mathcal {U}(\log (200), \log (10\,000))$$ in order to cover the allowed PBH halo parameter space. Figure 6. View largeDownload slide Same as Fig. 5, but for a Dehnen halo with a cusp (γ = 1) at t = 0. The prior distribution of log (R0, DM) is rescaled to $${\sim } \mathcal {U}(\log (200), \log (10\,000))$$ in order to cover the allowed PBH halo parameter space. Figure 7. View largeDownload slide Same as Fig. 6, but for a NFW halo truncated at 10 kpc at t = 0. Figure 7. View largeDownload slide Same as Fig. 6, but for a NFW halo truncated at 10 kpc at t = 0. Table 2. Results from 105 MCMC realizations assuming cored and cuspy DM density profiles. Model  log(R0, DM)  log(R0, *)  log(MDM)  log(MBH)  Denen core  $$2.77^{+0.17}_{-0.22}$$  $$1.47^{+0.05}_{-0.23}$$  $$9.12^{+0.42}_{-0.61}$$  $$0.93^{+0.23}_{-0.54}$$  Denen cusp  $$3.34^{+0.26}_{-0.28}$$  $$1.46^{+0.05}_{-0.27}$$  $$8.92^{+0.55}_{-0.44}$$  $$0.78^{+0.11}_{-0.41}$$  NFW cusp  $$3.15^{+0.45}_{-0.38}$$  $$1.46^{+0.043}_{-0.278}$$  $$9.23^{+0.33}_{-0.62}$$  $$0.76^{+0.14}_{-0.38}$$  Model  log(R0, DM)  log(R0, *)  log(MDM)  log(MBH)  Denen core  $$2.77^{+0.17}_{-0.22}$$  $$1.47^{+0.05}_{-0.23}$$  $$9.12^{+0.42}_{-0.61}$$  $$0.93^{+0.23}_{-0.54}$$  Denen cusp  $$3.34^{+0.26}_{-0.28}$$  $$1.46^{+0.05}_{-0.27}$$  $$8.92^{+0.55}_{-0.44}$$  $$0.78^{+0.11}_{-0.41}$$  NFW cusp  $$3.15^{+0.45}_{-0.38}$$  $$1.46^{+0.043}_{-0.278}$$  $$9.23^{+0.33}_{-0.62}$$  $$0.76^{+0.14}_{-0.38}$$  View Large A comparison of the constraints using the three different halo profiles shows that the results are very close to each other, and that these galaxies are DM-dominated: the posterior distribution of the central DM density lies in a narrow range of 1–2 M⊙pc − 3 while the stellar density is two orders of magnitude lower. We also experimented with models having no DM component at all. In this case, the velocity dispersion of stars would be much lower than the observed value for the given combination of total stellar mass and half-mass radius. By contrast, globular clusters have a far larger stellar mass while occupying the same region in the R*–σ plane. The existence of DM also ensures the logic of using ultra-faint dwarf galaxies to constrain PBH mass is self-consistent. The constraint on PBH mass in the massive end comes from a maximum heating rate, which is more stringent than that from the lower mass end, mainly due to the assumptions upon which the stellar component was initialized, as we discuss in Section 4. This can also be seen in the contour of R*, 0 and log (MBH), which is a continuous distribution. Models with R*, 0 between 20 and 30 pc are certainly allowed, with a modest heating from PBH that brings the stellar size and velocity dispersion close to the observed values. Note that the peak of MBH does not correspond to the peak of R*, 0. Therefore, the likelihood surface in the plane of R0, * and MBH is more informative, and the constraints on MBH is limited by our prior knowledge of R0, *. The estimate can be improved if we have better prior information about the initial size of the stellar component. There is a clear dichotomy in the spread in metallicity between dwarf galaxies and globular clusters of the same luminosities (see Willman & Strader 2012, most of the observed globular clusters show <0.1 dex metallicity spread). The large metallicity spread (>0.2 dex) among the member stars in the dwarf galaxies could impose more stringent constraints on the initial stellar component size. Gas confined in a small volume would lead to efficient mixing due to multiple supernova events, hence a more uniform metallicity distribution is incompatible with observations. If the initial size R0, * is above 20 pc, then MBH would fall below 10 M⊙. 4 DISCUSSIONS 4.1 Model limitations In our models, the strongest prior is the total DM halo mass, which fulfills the requirement that these haloes should be just slightly more massive than the minimum mass required for star formation. For the other three free parameters, we have adopted rather agnostic prior distributions. The initial distribution of stellar component, which follows a Plummer sphere, is a reasonable assumption. Our model assumes that DM halo has a density core. This is a reasonable assumption given that the collisional effect on a cusped density profile is well understood. To demonstrate this, we have run two additional sets of MCMC realizations with a Dehnen density cusp (γ = 1) at t = 0 and an NFW cuspy profile. As listed in Table 2 and illustrated in Figs 5–7, the final distributions of the four parameters are very close among all three models, although the peak of log (MBH) drops slightly from 0.93 to 0.76. This is due to extra heating as the central velocity dispersion is low for a cuspy density profile. Our choice of prior distribution of total DM mass is also a reasonable one. As illustrated in Fig. 5, the total DM halo mass MDM does not have strong correlation with either MBH or R*, 0. The covariance between MDM and R0, DM covers a region defined by a constant DM density in the range of 1–2 M⊙pc − 3. In contrast, PBH mass is strongly correlated with the initial size of the stellar component. Therefore, our constraints on PBH mass is unaffected by the prior distribution of total DM mass. Moreover, if more compact ultra-faint dwarfs are to be discovered in the near future, their cumulative number could impose yet another constraint against the possibility of massive M⊙ PBHs as DM. There are several limitations in our approach that need to be stressed. First, we treat the haloes in isolation while in reality they are in highly tumultuous galactic environments, subjected to various mass loss processes (Zhu et al. 2016). Secondly, we use the measurements of the compact ultra-faint dwarfs as an input for the likelihood function and treat the reported uncertainties as the widths of Gaussians in the likelihood function. We have used a smaller half-mass size for Ret II and Hor I between those reported by Koposov et al. (2015) and Bechtol et al. (2015). Accurate determination of the membership of stars in such faint galaxies is an extremely challenging task. In addition, unresolved binary stars may significantly inflate the observed velocity dispersion in these cold systems, although they are unlikely to completely dominate the measurements (McConnachie & Côté 2010). Should the true size of these stellar systems change, our modelling and conclusions need to be revisited accordingly. In our model, there is a built-in bias which puts more emphasis on the more massive PBH for the following reason. If the collisional relaxation effect is negligible, only a limited region around rh = 40 pc and σ = 4 km s−1 is allowed. On the other hand, if the model starts with a smaller initial size (with a smaller velocity dispersion), a modest heating will bring both the size and the velocity dispersion closer to the observed values. Therefore, some heating effect can be slightly more favoured than no heating effect if there are more models with modest heating allowed than those with little heating. 4.2 Comparison with previous works Our analysis improves upon the previous calculations by Brandt (2016) and Koushiappas & Loeb (2017) by using a FP code and an MCMC approach. In particular, the FP code self-consistently evolved the two components system with two mass species. In the previous studies, the velocity dispersion profiles of the two components are assumed to be the same, while our modelling with the FP code clearly shows this is not the case. Moreover, our calculations show that the heating rate from the massive species on the light species is slightly less efficient than t0.5, a result of partial equipartition of energy due to the dominant contribution from PBHs. In addition, the combination of the FP code and an MCMC approach is able to utilize both size and kinematics of the observed ultra-compact dwarf galaxies by efficiently sampling the parameter space. We note that the stellar density profile by the end of the integration contains a ‘diffused’ core in our FP calculation. This is not the same as the ring profile predicted by Koushiappas & Loeb (2017). In Fig. 8, we demonstrate the effect of PBH mass on the evolution of stellar surface density. The stars initially distribute in a Plummer sphere with a mass of M* = 350 M⊙ and a scale radius of R0, * = 20 pc. The evolved projected surface density profiles of stars at t = 12 Gyr are shown for different PBH mass ranges from MBH = 3–100 M⊙. As the PBH mass increases, the two-body heating intensifies, which leads to a stronger evolution in the core, featuring a larger core size and a lower surface density in the central region while an increase in the outer part. Clearly, the size evolution in our simulations is global, while Koushiappas & Loeb (2017) assumes that each shell evolved independently. In addition, we use an isotropic velocity distribution, which is different from the highly circular orbits, an implicit assumption made in Koushiappas & Loeb (2017). These different assumptions and treatments result in different stellar surface density profiles between our simulations and the previous analytical calculations. Figure 8. View largeDownload slide Projected stellar surface density profile from the simulation at t = 12 Gyr with different PBH masses, as represented by different colours. The initial stellar population is set up as a Plummer sphere with M* = 350 M⊙ and R0, * = 20 pc. For $$\rm {M_{BH} = 3\, \rm {M_{\rm {\odot }}}}$$, there is little evolution of the density profile for stars. As the PBH mass increases, there is a stronger evolution in the core, resulting in a larger core size and a lower surface density in the central region. Figure 8. View largeDownload slide Projected stellar surface density profile from the simulation at t = 12 Gyr with different PBH masses, as represented by different colours. The initial stellar population is set up as a Plummer sphere with M* = 350 M⊙ and R0, * = 20 pc. For $$\rm {M_{BH} = 3\, \rm {M_{\rm {\odot }}}}$$, there is little evolution of the density profile for stars. As the PBH mass increases, there is a stronger evolution in the core, resulting in a larger core size and a lower surface density in the central region. Our model assumes that the DM is completely composed of PBHs. It was suggested by Brandt (2016) that PBH above 10 M⊙ may be firmly excluded as the primary candidate of DM in dwarf galaxies. The constraint from Koushiappas & Loeb (2017) is slightly stronger, with masses greater than 6 M⊙ excluded. The constraint derived from our analysis is slightly less restrictive compared to these two studies. 4.3 Observational signatures of PBHs with millisecond time delay Our analysis is consistent with the suggestion that PBHs in the mass range of 25–100 M⊙ are ruled out as the main component of DM (Brandt 2016; Carr et al. 2017). If the LIGO detections are from mergers of PBH binaries of ∼ 30 M⊙, it is only possible if they are in the massive tail of an extended PBH mass function (Magee & Hanna 2017). However, any PBHs with a mass substantially higher than 10 M⊙ would violate the constraints from the dynamical analysis of dwarf galaxies. In the near future, aLIGO observations may soon provide constraints on PBHs around 10 M⊙. It remains possible that some fraction of the DM may be PBHs in this mass range. Apart from the contribution to the total DM budget, the existence of PBHs will provide important constraints on certain inflation theories. In addition, gravitational lensing due to ∼10 M⊙ PBHs will induce a time delay on the order of milliseconds (Mao 1992; Muñoz et al. 2016). There are at least two known astrophysical phenomena that can produce detectable signals on this rapid time-scale. One is fast radio burst (Lorimer et al. 2007). If fast radio bursts have indeed cosmological origin (Spitler et al. 2016), then one would expect repeated bursts with millisecond delay from the same location of the sky within a large number (Fialkov & Loeb 2017) of bursts due to the PBH lenses (Muñoz et al. 2016). Another promising source in the kHz regime for the detection of PBHs would be gravitational waves from stellar mass BH mergers. Similar to the mechanism of femtolensing of γ-ray burst (Gould 1992), millisecond delays could induce a characteristic interference pattern in the detected waveforms in the merger phase due to the interference of two lensed images. Perhaps gravitational waves are better sources than radio bursts since they do not suffer from interstellar dispersions. This technique could be feasible with the next generation detectors such as the ‘Cosmic Explorer’ (Abbott et al. 2017). With a large number of gravitational wave detections, we could place a strong constraint on the mass fraction of PBHs in ΩM, since the lensing optical depth is at the same order as the mass fraction of PBHs for cosmological sources (Press & Gunn 1973). 5 SUMMARY In this study, we improve the recent analytical work of Brandt (2016) and Koushiappas & Loeb (2017) on the investigation of the possibility of DM consisting entirely of PBHs, by combining Fokker–Planck simulations with MCMC realizations. The former method, implemented in the code phaeflow, accurately follows the joint evolution of stars and PBH-DM driven by two-body relaxation, for a given choice of parameters. By combining it with the MCMC approach, we explore the parameter space and use the half-light radius and central stellar velocity dispersions of observed compact ultra-faint dwarf galaxies as joint priors in the MCMC sampling to constrain the PBH-DM parameters. Our findings can be summarized as follows: For a DM halo in the mass range of $$10^6\text{--}10^9\, \rm {M_{\rm {\odot }}}$$, if the DM is composed of PBHs entirely, collisional relaxation quickly transforms a cusp to a cored density profile. For a system with stars and PBH-DM, two-body relaxation between the stars and PBH-DM drives up the stellar core size, and heating from the PBHs increases the central stellar velocity dispersion. The size evolution of the stellar core has a growth rate of rh ∝ t0.4, which is slightly slower than that from the full energy equipartition (rh ∝ t0.5). The joint constraints from the observed half-light radius and central stellar velocity dispersions of compact ultra-faint dwarf galaxies leave a narrow parameter space of DM haloes consist of PBHs. Results from the MCMC simulations suggest that the central density of DM halo has a narrow range of 1–2 M⊙pc − 3, and the mass of PBHs is most likely within 2–14 M⊙, if the DM is completely made of PBH. We note that our constraint of PBHs is slightly less restrictive than those of Brandt (2016) and Koushiappas & Loeb (2017). Although PBH with mass above 10 M⊙ can be excluded as the primary candidate of DM from our results, we emphasize that our models are inconclusive, and that the nature of DM remains elusive. ACKNOWLEDGEMENTS We thank the anonymous referee for a constructive report which has helped improve our paper. It is our great pleasure to thank Avi Loeb, Lars Hernquist, Savvas Koushiappas and Peter Mészáros for stimulating discussions. YL acknowledges support from NSF grants AST-0965694, AST-1009867, AST-1412719, and MRI-1626251. The numerical computations and data analysis in this paper have been carried out on the CyberLAMP cluster supported by MRI-1626251, operated and maintained by the Institute for CyberScience at the Pennsylvania State University, as well as the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. The Institute for Gravitation and the Cosmos is supported by the Eberly College of Science and the Office of the Senior Vice President for Research at the Pennsylvania State University. EV acknowledges support from the European Research Council under the 7th Framework Programme (grant 321067). Footnotes 1 phaseflow is part of the AGAMA library for galaxy modelling, available from https://github.com/GalacticDynamics-Oxford/Agama/. REFERENCES Abbott B. P. et al.  , 2017, Class. Quantum Gravity , 34, 044001 https://doi.org/10.1088/1361-6382/aa51f4 CrossRef Search ADS   Afshordi N., McDonald P., Spergel D. N., 2003, ApJ , 594, L71 https://doi.org/10.1086/378763 CrossRef Search ADS   Amorisco N. C., 2017, AJ , 844, 16 CrossRef Search ADS   Bechtol K. et al.  , 2015, ApJ , 807, 50 https://doi.org/10.1088/0004-637X/807/1/50 CrossRef Search ADS   Bianchini P., van de Ven G., Norris M. A., Schinnerer E., Varri A. L., 2016, MNRAS , 458, 3644 https://doi.org/10.1093/mnras/stw552 CrossRef Search ADS   Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition . Princeton Univ. Press, Princeton, NJ Bird S., Cholis I., Muñoz J. B., Ali-Haïmoud Y., Kamionkowski M., Kovetz E. D., Raccanelli A., Riess A. G., 2016, Phys. Rev. Lett. , 116, 201301 https://doi.org/10.1103/PhysRevLett.116.201301 CrossRef Search ADS PubMed  Brandt T. D., 2016, ApJ , 824, L31 https://doi.org/10.3847/2041-8205/824/2/L31 CrossRef Search ADS   Carr B. J., Hawking S. W., 1974, MNRAS , 168, 399 https://doi.org/10.1093/mnras/168.2.399 CrossRef Search ADS   Carr B., Kühnel F., Sandstad M., 2016, Phys. Rev. D , 94, 083504 https://doi.org/10.1103/PhysRevD.94.083504 CrossRef Search ADS   Carr B., Raidal M., Tenkanen T., Vaskonen V., Veermäe H., 2017, Phys. Rev. D , 96, 023514 https://doi.org/10.1103/PhysRevD.96.023514 CrossRef Search ADS   Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS , 454, 2981 https://doi.org/10.1093/mnras/stv2165 CrossRef Search ADS   Contenta F. et al.  , 2017, MNRAS , preprint (arXiv:1705.01820) Dehnen W., 1993, MNRAS , 265, 250 https://doi.org/10.1093/mnras/265.1.250 CrossRef Search ADS   Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Stinson G. S., Knebe A., 2014, MNRAS , 441, 2986 https://doi.org/10.1093/mnras/stu729 CrossRef Search ADS   Diemer B., Kravtsov A. V., 2015, ApJ , 799, 108 https://doi.org/10.1088/0004-637X/799/1/108 CrossRef Search ADS   Fialkov A., Loeb A., 2017, ApJ , 846, L27 https://doi.org/10.3847/2041-8213/aa8905 CrossRef Search ADS   Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP , 125, 306 https://doi.org/10.1086/670067 CrossRef Search ADS   Frampton P. H., Kawasaki M., Takahashi F., Yanagida T. T., 2010, J. Cosmol. Astropart. Phys. , 4, 023 https://doi.org/10.1088/1475-7516/2010/04/023 CrossRef Search ADS   Garrison-Kimmel S., Bullock J. S., Boylan-Kolchin M., Bardwell E., 2017, MNRAS , 464, 3108 https://doi.org/10.1093/mnras/stw2564 CrossRef Search ADS   Gnedin N. Y., Kaurov A. A., 2014, ApJ , 793, 30 https://doi.org/10.1088/0004-637X/793/1/30 CrossRef Search ADS   Gould A., 1992, ApJ , 386, L5 https://doi.org/10.1086/186279 CrossRef Search ADS   Hawking S., 1971, MNRAS , 152, 75 https://doi.org/10.1093/mnras/152.1.75 CrossRef Search ADS   Hernquist L., 1990, ApJ , 356, 359 https://doi.org/10.1086/168845 CrossRef Search ADS   Hinton S., 2016, JOSS , 1, 45 CrossRef Search ADS   Hirano S., Sullivan J. M., Bromm V., 2017, MNRAS , 473, L6 CrossRef Search ADS   Hu W., Barkana R., Gruzinov A., 2000, Phys. Rev. Lett. , 85, 1158 https://doi.org/10.1103/PhysRevLett.85.1158 CrossRef Search ADS PubMed  Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, Phys. Rev. D , 95, 043541 https://doi.org/10.1103/PhysRevD.95.043541 CrossRef Search ADS   Kashlinsky A., 2016, ApJ , 823, L25 https://doi.org/10.3847/2041-8205/823/2/L25 CrossRef Search ADS   Khlopov M. Y., 2010, Res. Astron. Astrophys. , 10, 495 https://doi.org/10.1088/1674-4527/10/6/001 CrossRef Search ADS   Koposov S. E., Belokurov V., Torrealba G., Evans N. W., 2015, ApJ , 805, 130 https://doi.org/10.1088/0004-637X/805/2/130 CrossRef Search ADS   Koushiappas S. M., Loeb A., 2017, Phys. Rev. Lett. , 119, 041102 https://doi.org/10.1103/PhysRevLett.119.041102 CrossRef Search ADS PubMed  Lorimer D. R., Bailes M., McLaughlin M. A., Narkevic D. J., Crawford F., 2007, Science , 318, 777 https://doi.org/10.1126/science.1147532 CrossRef Search ADS PubMed  Ma X. et al.  , 2017, MNRAS , preprint (arXiv:1706.06605) Magee R., Hanna C., 2017, AJ , 845, 5 CrossRef Search ADS   Mao S., 1992, ApJ , 389, L41 https://doi.org/10.1086/186344 CrossRef Search ADS   McConnachie A. W., Côté P., 2010, ApJ , 722, L209 https://doi.org/10.1088/2041-8205/722/2/L209 CrossRef Search ADS   Meszaros P., 1975, A&A , 38, 5 Muñoz J. B., Kovetz E. D., Dai L., Kamionkowski M., 2016, Phys. Rev. Lett. , 117, 091301 https://doi.org/10.1103/PhysRevLett.117.091301 CrossRef Search ADS PubMed  Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ , 490, 493 https://doi.org/10.1086/304888 CrossRef Search ADS   Okamoto T., Gao L., Theuns T., 2008, MNRAS , 390, 920 https://doi.org/10.1111/j.1365-2966.2008.13830.x CrossRef Search ADS   Plummer H. C., 1911, MNRAS , 71, 460 https://doi.org/10.1093/mnras/71.5.460 CrossRef Search ADS   Pontzen A., Governato F., 2014, Nature , 506, 171 https://doi.org/10.1038/nature12953 CrossRef Search ADS PubMed  Press W. H., Gunn J. E., 1973, ApJ , 185, 397 https://doi.org/10.1086/152430 CrossRef Search ADS   Quinlan G. D., 1996, New A , 1, 255 https://doi.org/10.1016/S1384-1076(96)00018-8 CrossRef Search ADS   Read J. I., Agertz O., Collins M. L. M., 2016, MNRAS , 459, 2573 https://doi.org/10.1093/mnras/stw713 CrossRef Search ADS   Sawala T. et al.  , 2015, MNRAS , 448, 2941 https://doi.org/10.1093/mnras/stu2753 CrossRef Search ADS   Shapiro P. R., Iliev I. T., Raga A. C., 1999, MNRAS , 307, 203 https://doi.org/10.1046/j.1365-8711.1999.02609.x CrossRef Search ADS   Spitler L. G. et al.  , 2016, Nature , 531, 202 https://doi.org/10.1038/nature17168 CrossRef Search ADS PubMed  Spitzer L., Jr, 1969, ApJ , 158, L139 https://doi.org/10.1086/180451 CrossRef Search ADS   The LIGO Scientific Collaboration, 2016a, Phys. Rev. Lett. , 116, 061102 https://doi.org/10.1103/PhysRevLett.116.061102 CrossRef Search ADS   The LIGO Scientific Collaboration, 2016b, Phys. Rev. Lett. , 116, 241103 https://doi.org/10.1103/PhysRevLett.116.241103 CrossRef Search ADS   The LIGO Scientific Collaboration, 2017a, Phys. Rev. Lett. , 119, 141101 CrossRef Search ADS   The LIGO Scientific Collaboration, 2017b, Phys. Rev. Lett. , 118, 221101 https://doi.org/10.1103/PhysRevLett.118.221101 CrossRef Search ADS   Trenti M., van der Marel R., 2013, MNRAS , 435, 3272 https://doi.org/10.1093/mnras/stt1521 CrossRef Search ADS   Vasiliev E., 2017, ApJ , 848, 10 https://doi.org/10.3847/1538-4357/aa8cc8 CrossRef Search ADS   Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS , 423, 3740 https://doi.org/10.1111/j.1365-2966.2012.21182.x CrossRef Search ADS   Vogelsberger M., Zavala J., Cyr-Racine F.-Y., Pfrommer C., Bringmann T., Sigurdson K., 2016, MNRAS , 460, 1399 https://doi.org/10.1093/mnras/stw1076 CrossRef Search ADS   Weinberg D. H., Bullock J. S., Governato F., Kuzio de Naray R., Peter A. H. G., 2015, Proc. Natl. Acad. Sci. , 112, 12249 https://doi.org/10.1073/pnas.1308716112 CrossRef Search ADS   Wheeler C., Oñorbe J., Bullock J. S., Boylan-Kolchin M., Elbert O. D., Garrison-Kimmel S., Hopkins P. F., Kereš D., 2015, MNRAS , 453, 1305 https://doi.org/10.1093/mnras/stv1691 CrossRef Search ADS   Willman B., Strader J., 2012, AJ , 144, 76 https://doi.org/10.1088/0004-6256/144/3/76 CrossRef Search ADS   Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., Muñoz R. R., Simon J. D., Avedo F. F., 2010, MNRAS , 406, 1220 Yoshida N., Sokasian A., Hernquist L., Springel V., 2003, ApJ , 591, L1 https://doi.org/10.1086/376963 CrossRef Search ADS   Zhao H., 1996, MNRAS , 278, 488 https://doi.org/10.1093/mnras/278.2.488 CrossRef Search ADS   Zhu Q., Marinacci F., Maji M., Li Y., Springel V., Hernquist L., 2016, MNRAS , 458, 1559 https://doi.org/10.1093/mnras/stw374 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations