Improved constraints from ultra-faint dwarf galaxies on primordial black holes as dark matter

Soon after the recent first ever detection of gravitational waves from merging black holes it has been suggested that their origin is primordial. Appealingly, a sufficient number of primordial black holes (PBHs) could also partially or entirely constitute the dark matter (DM) in the Universe. However, recent studies on PBHs in ultra-faint dwarf galaxies (UFDGs) suggest that they would dynamically heat up the stellar component due to two-body relaxation processes. From the comparison with the observed stellar velocity dispersions and the stellar half-light radii it was claimed that only PBHs with masses $\lesssim10\,M_\odot$ can significantly contribute to the DM. In this work, we improve the latter constraints by considering the largest observational sample of UFDGs and by allowing the PBH masses to follow an extended (log-normal) distribution. By means of collisional Fokker-Planck simulations, we explore a wide parameter space of UFDGs containing PBHs. The analysis of the half-light radii and velocity dispersions resulting from the simulations leads to three general findings that exclude PBHs with masses $\sim\mathcal{O}(1$-$100)\,M_\odot$ from constituting all of the DM: (i) We identify a critical sub-sample of UFDGs that only allows for $\sim\mathcal{O}(1)\,M_\odot$ PBH masses; (ii) for any PBH mass, there is an UFDG in our sample that disfavours it; (iii) the spatial extensions of a majority of simulated UFDGs containing PBHs are too large to match the observed.


INTRODUCTION
Compelling evidence for the existence of dark matter (DM) has been found on almost all astrophysical scales; still, its very nature remains elusive.Solving the DM puzzle has therefore become one of the greatest challenges in presentday astrophysics.Recently, the first direct detections of gravitational waves released in black hole binary mergers (Abbott et al. 2019) revived interest in the conjecture that primordial black holes (PBHs) may constitute partly or entirely the DM in the Universe (e.g.Clesse & García-Bellido 2018;Bird et al. 2016;Sasaki et al. 2016;Carr et al. 2016;Raidal et al. 2017;García-Bellido 2017;Clesse & García-Bellido 2017).That being the case, numerous PBHs would have formed out of sufficiently large overdensities in the primordial matter power spectrum; such overdensities would eventually have collapsed during the radiation-dominated E-mail: StegmannJ@cardiff.ac.uk epoch, producing a population of PBHs that survived to the present (Rice & Zhang 2017).The PBH masses at time of formation are expected to be roughly the enclosed mass within the horizon at that time (Sasaki et al. 2018).For this reason, depending on the precise formation time, PBHs can span a wide range of masses.Due to Hawking (1974) radiation, PBHs with masses M ∼ < 10 −19 M should have already evaporated by now (Page 1976), whereas heavier PBHs could in principle have survived until today, contributing to the present-day DM.Amongst those, PBHs with masses O(1-100) M are of particular interest as they can account for the astrophysical origin of BH binaries detected through gravitational waves.
For these PBHs, several observations constrain the possibility that they constitute a significant fraction of the DM.Firstly, if a significant fraction of the Galactic DM halo consists of PBHs, the light coming from background sources would occasionally exhibit a microlensing pattern (Paczynski 1986).That is, PBHs that intersect the optical axis be-tween observer and source deflect its light rays, resulting in a characteristic magnification of the apparent luminosity over a finite period of time.Several missions aimed at looking for such microlensing events on nearby stars and distant quasars have been carried out in the last decades (e.g.Niikura et al. 2019;Tisserand et al. 2007;Alcock et al. 2001;Mediavilla et al. 2017).In total, a few dozen microlensing events have been observed that last up to several months and whose rates put stringent upper limits on PBHs with masses up to O(10) M as a dominant DM component (Carr et al. 2017).Secondly, a PBH would accrete baryons getting sufficiently close to its horizon; if this happens in the early Universe, the resulting radiation would leave signatures on the power spectrum of the cosmic microwave background (CMB) and on its polarisation anisotropies (Ali-Haïmoud & Kamionkowski 2017; Sasaki et al. 2018).From this, conservative estimates suggest PBHs with masses roughly above O(100) M to be excluded from contributing significantly to the DM (Carr et al. 2017); however, one should keep in mind that such constraints intrinsically suffer from uncertainties in the modelling of the accretion physics.In this paper, we investigate whether observations of ultrafaint dwarf galaxies (UFDGs) are compatible with the O(1-100) M PBH scenario.Owing to their extraordinarily high ratio between dynamical mass and stellar mass, UFDGs represent an excellent candidate to test DM models (Simon 2019).Previous studies of PBHs in UFDGs typically investigated the dynamical imprint of PBHs on the stellar population.Brandt (2016) showed that a recently discovered star cluster in Eridanus II as well as the entire stellar populations of dwarf galaxies are vulnerable to dynamical heating by PBHs; based on the survival of the star cluster and that of a selected sub-sample of UFDGs, and assuming a monochromatic PBH mass function (i.e. a mass function in which all PBHs have the same mass), he put stringent upper bounds on PBHs with masses respectively ∼ > 5 and ∼ > 10 M as constituting the bulk of DM. Green (2016) generalised these results to an extended (log-normal) PBH mass function by studying the consequences on the Eridanus II star cluster and on the UFDGs' stellar populations.From this, she excluded broad PBH mass functions with a significant amount of O(10) M PBHs.In a different study, Koushiappas & Loeb (2017) investigated effects of mass segregation between PBHs and stars in Segue I.They specifically predicted a depletion of stars in the UFDG's centre as well as a ring structure in the projected stellar surface density profile.By comparing these predictions with the observation of Segue I, they excluded PBHs of single mass O(10) M from constituting a significant DM fraction.Finally, based on a similar sub-sample of UFDGs used by Brandt (2016), Zhu et al. (2018) subsequently aimed to reconstruct the UFDGs' stellar observables via a large number of galaxy simulations with varying the (monochromatic) PBH mass.From the (dis)agreement of observations with the related simulations, they concluded that only PBHs with masses 2-14 M can significantly contribute to the DM.Their simulations also appeared in conflict with the predictions from Koushiappas & Loeb (2017), as they did not exhibit any ring structure in the stellar density profile.In summary, dynamical constraints from UFDGs severely challenge the possibility that PBHs would constitute a sig-nificant fraction of the DM.Nevertheless, there is room left for improvement of the previous results that could strengthen such dynamical constraints.For instance, previous studies are typically restricted to a limited sample of UFDGs.In fact, only the very faint end of UFDGs has been considered which allows PBH masses around O(1) M .In contrast, the growing number of detections invites us to take a larger sample into account.A priori, one would expect three possible outcomes for this more complete investigation.Either the added UFDGs allow the same PBH masses, they allow a different mass window, or they exclude PBHs of any mass from constituting a significant part of the DM.Evidently, the latter two results would severely challenge the PBH scenario.On the other hand, with the exception of Green (2016), previous studies just considered monochromatic PBH mass functions.However, as pointed out for instance by Carr et al. (2017) and Clesse & García-Bellido (2015), this assumption is neither realistic with respect to most of the proposed PBH formation mechanisms nor does it usually lead to the same constraints.Hence, the aim of this work is to overcome this insufficiency by providing a general constraint on the existence of PBHs in UFDGs.For this purpose, we follow the approach of Zhu et al. (2018) but we consider a most extensive sample of UFDGs and an extended PBH mass function.Specifically, we investigate whether PBHs can constitute all of the DM (hereafter, PBH-DM).This paper is organised as follows: Section 2 is dedicated to the description of the methods we use.In Section 3, we present the results that are discussed in Section 4.

METHODS
The primary objective of this work is to investigate whether the stellar properties of observed UFDGs are compatible with DM entirely constituted by PBHs.Our analysis is motivated by the fact that the relaxation time can become comparable to the Hubble time in this scenario.Consequently, PBH-DM in UFDGs has to be treated as a collisional fluid potentially affecting the shape of the UFDGs.Such collisional evolution will be reflected in the observable stellar populations.For this purpose, we perform large sets of 2 × 10 5 collisional simulations for each UFDG in our sample (described in Section 2.1) and we compare the resulting observational properties with the observed ones.Each simulation models the interplay between stars and PBH-DM altogether specified by a set of input parameters.Furthermore, as opposed to previous work, we do not restrict our analysis to a monochromatic PBH mass function, but also investigate the consequences of an extended (log-normal) PBH mass function (see Section 2.2) by explicitly including its model parameters as input parameters.From comparison between simulations and observations, we are thus able to infer the probabilities for each input parameter set to reproduce the stellar properties, i.e. we evaluate if a given set of parameters is reproduced in the associated UFDG or not.Finally, the combined analysis for all UFDGs allows us to impose a generic constraint on PBH-DM.More specifically, given an UFDG, let D = {dj} and Θ = {θi} be the set of observed stellar properties and the set of Improved constraints from UFDGs on PBHs as DM 3 input parameters for one simulation, respectively.Then we draw 2 × 10 5 samples from the posterior distribution, where p(Θ) and p(D | Θ) are the prior distribution and the likelihood function, respectively, by the use of emcee  (Bahcall & Wolf 1976;Vasiliev 2017), and monochromatic PBH mass functions (Zhu et al. 2018).
In the following, we explicate each aspect relevant for our approach, including the observational sample of UFDGs (Section 2.1), the selection of the PBH mass functions (Section 2.2), the models of UFDGs (Section 2.3), and the choice of the prior distribution p(Θ) as it appears in Eq. (1) (Section 2.4).

Sample of observed UFDGs
Following the definition of Simon (2019), UFDGs are dwarf galaxies with V-band luminosities LV ≤ 10 5 L , corresponding to absolute magnitudes MV ≥ −7.7.Considering their stellar masses, surface brightnesses, sizes, dynamical masses, and metallicities, the hitherto observed UFDGs do not appear as fundamentally distinct from other dwarf galaxies (Simon 2019).Unlike in previous studies on PBH-DM in UFDGs, we therefore do not restrict ourselves to a selected, seemingly distinct sub-sample of UFDGs (as done in, e.g.Brandt 2016;Zhu et al. 2018), but instead consider an extensive sample of 27 UFDGs whose projected two-dimensional stellar half-light radii (r h, ) and line-of-sight stellar velocity dispersions (σ ) have been observationally constrained.We list them in Table 1.We assume that each star in the simulation shares the same mass and the same mass-to-light ratio (cf.Section 2.3).Consequently, the (projected twodimensional) half-light radius is equivalent to the (projected two-dimensional) half-mass radius of the stars.In order to resolve both the stellar kinematics and their structural properties, we base the comparison between simulations and observations on the stellar velocity dispersion and the deprojected stellar half-mass radius R h, = (4/3)r h, (Wolf et al. 2010), i.e.D = {σ , R h, }.Thus, if a simulation of a given UFDG observed with (σ ,D) Table 1) outputs the stellar velocity dispersion σ ,S and the deprojected half-mass radius R h, ,S , its likelihood function as it appears in Eq. ( 1) is evaluated to where SN denotes the split-normal distribution.The σi (i = 1-4) are the respective measurement uncertainties.This approach accounts for the asymmetry in some of the reported measurements.For those UFDGs where only an upper limit for σ is given, SN in Eq. ( 2) is replaced by a half-normal distribution whose parameter is chosen such that its cumulative distribution function evaluated at the given upper limit matches the respective confidence level.Lastly, when analysed, the UFDGs' stellar populations turn out to be fairly old (t ≈ 9-14 Gyr).Therefore, we run each simulation over t = 12 Gyr.

Extended PBH mass functions
The two key improvements of this paper compared to Zhu et al. (2018) are the generalisations to a severely enlarged sample of UFDGs, presented in Section 2.1, and to PBHs covering an extended mass range.Usually, the latter are described by means of a mass function ψ(M ), which determines their mass density ρPBH and number density nPBH: Here, ρDM denotes the total DM mass density.So defined, the mass function is normalised to the total fraction ftot ≡ ρPBH/ρDM: Evidently, 0 ≤ ftot ≤ 1 must hold.In this paper, we investigate whether PBHs can constitute all of the DM, i.e. whether ftot = 1 or not.1 Henceforth, we set ftot = 1 and ρPBH and ρDM can be used equivalently.By mean of Eqs (3), (4), and (5), we introduce the mass function ψ(M ) such that it agrees with the notation of Carr et al. (2017), who comprehensively generalised previously existing constraints on monochromatic mass functions to extended.In the following, we list the functional forms of both the monochromatic mass function and the log-normal mass function, which we investigate as a commonly discussed representative for generic extended mass functions.

UFDG Name
Improved constraints from UFDGs on PBHs as DM 5 • A monochromatic mass function: where δ denotes the Dirac delta function.By definition, this mass function only depends on the PBH mass M .Typical mass functions emerging from a critical collapse of density fluctuations with a δ-function power spectrum are relatively narrow and can therefore also be well approximated with a monochromatic mass function (Carr et al. 2017).When performing the simulations, we adopt the PBH mass as a free input parameter Mc ∈ Θ.
• A log-normal mass function: where µ > 0 and σ > 0. The mean PBH mass in the assumption of a log-normal mass function is M = µ exp −σ 2 /2 .Such mass function is a good approximation for a large class of PBH formation scenarios, e.g.axion-curvaton, running-mass, and single field double inflation (Green 2016;Kannike et al. 2017).The (µ, σ) pairs allowed in the window relevant for the operational gravitational wave detectors [O(1-100) M ; Abbott et al. 2019] are mainly constrained by (i) the results from searches for microlensing events on stars in our Galactic neighbourhood and (ii) the effect PBH gas accretion would have on the CMB temperature and ionization history (Sasaki et al. 2018;Carr et al. 2017).Together, the allowed points roughly constitute a sub-plane described as (see figure 3 in Carr et al. 2017): It is worth mentioning that these points are ruled out completely if one further takes the following constraints into account (Carr et al. 2017): (iii) the survival of the stellar cluster in Eridanus II and of the entire stellar populations in UFDGs (Brandt 2016;Koushiappas & Loeb 2017;Green 2016;Zhu et al. 2018) and (iv) the survival of wide binaries in the Milky Way (Monroy-Rodríguez & Allen 2014).These constraints are somewhat less restrictive as they rely on further astrophysical assumptions (Carr et al. 2017).For instance, the stellar cluster in Eridanus II could have only recently spiralled down into the centre of the galaxy, where dynamical heating by PBHs becomes effective (Brandt 2016).Additionally, wide binaries are in principle hard to detect, yielding some uncertainty in identifying them and consequently in drawing any conclusion on the PBH abundance (Sasaki et al. 2018).Moreover, we emphasise that, most recently, even previous microlensing constraints were called into question as spatial PBH clustering (García-Bellido & Clesse 2018) and updated galactic rotation curves (Hawkins 2015) tend to relax them.When we implement the log-normal mass function in the simulation, we take both parameters, µ and σ, as free input parameters: µ, σ ∈ Θ.

Initial UFDG configurations
Each point in the MCMC represents an UFDG that consists of stars and PBHs and whose evolution is simulated over an integration time t = 12 Gyr with PhaseFlow.The numerical value for t is motivated in Section 2.1 by the given observational UFDG sample.For any simulation, we assume that the initial radial density distribution of the PBHs is described by a Dehnen (1993) sphere: where MDM, R0,DM, and γ are the DM total mass, scale radius, and asymptotic inner slope, respectively.Whilst taking MDM and R0,DM as free input parameters (MDM, R0,DM ∈ Θ), we fix γ = 0 for all simulations, as PBHs have been shown (Zhu et al. 2018) to be a cusp-core transformer by rapidly converting an initially cuspy (γ = 1) density profile into a cored one (γ = 0).We identified the same behaviour in direct N -body simulations of dwarf galaxies with PBH-DM that is initially described by a cuspy density profile.Most recently, this result was also found by Boldrini et al. (2019).In their simulations, Zhu et al. (2018) implement a single DM component consisting of monochromatic PBHs that follow Eq. ( 9).Here we also investigate the extended mass functions introduced in Section 2.2 by approximating them with N = 50 mass bins, each of which corresponds to one component in the simulations.We emphasise that this generalisation is a main improvement presented in this paper since it allows us to constrain physically more realistic PBH mass functions than the monochromatic one.More specifically, given a mass function ψ(M ), is the cumulative DM fraction of PBHs with masses M ∈ [0, M th ] [cf.with Eq. ( 5)].For all mass functions introduced in Section 2.2, F (0; M th ) is invertible so that there are uniquely determinable mass thresholds M th,0 and M th,1 at which F (0; M th,0 ) = 0.5 per cent and F (0; M th,1 ) = 99.5 per cent, respectively.Once M th,0 and M th,1 are calculated, we divide the mass range in log 10space, [log 10 M th,0 , log 10 M th,1 ], in N bins of same length and assume that all PBHs in one bin have a mass equal to the bin's average, where M bin,0 and M bin,1 are the bin's upper and lower boundary, respectively.Consequently, the total DM mass contained within a bin is computed as its fractional mass, MDMF (M bin,0 ; M bin,1 ).Loosely speaking, we thus perform the transition from a monochromatic to an extended mass function by replacing the single monochromatic PBH-DM component with N PBH-DM components, each of which is in turn monochromatic but at different masses dictated by the given mass function.In this approach, we artificially neglect 1 per cent at the tails of the mass function, primarily in order to avoid bins with an average PBH mass greater than the bin's total DM mass, which would otherwise lead to numerical errors in the simulation.In Figure 1, this mass binning scheme is visualised for an exemplary mass function.
Concerning the initial distribution of the stars, we set up a Plummer (1911) sphere, described by where M and R0, are the total stellar mass and scale radius, respectively.Here, R0, is taken as a free input parameter (R0, ∈ Θ).The total stellar mass M is estimated from the luminosity of each UFDG separately (cf.Table 1) by assuming a stellar mass-to-light ratio of 2 M /L , since the UFDGs' stars are fairly old (Simon 2019).We also perform comparison runs with 1 and 3 M /L .Additionally, the stellar mass distribution is assumed to be monochromatic at 1 M .Even though assuming a single mass value for all the stars is a simplification, we do not expect it to impact our results as (i) stellar masses cover a less broad mass range compared to the one tested for PBHs and (ii) the low mass of stars, compared to the PBH mass obtained in this study, implies that the overall dynamical evolution (relaxation) of the simulated UFDGs is virtually not influenced by the stellar population.Finally, the average velocity changes per unit time experienced by the stars and PBHs due to mutual gravitational interactions are typically described by means of diffusion coefficients.These scale with the Coulomb logarithm ln Λ, where Λ = bmax/b90 is the ratio between the maximum impact parameter, bmax, and the impact parameter, b90, for which 90 degree scattering occurs (Binney & Tremaine 2008).Unfortunately, there is no precise definition for either b90 or bmax, as they can be somehow position-dependent for a given system.In fact, ln Λ is not straightforwardly defined even for the simplest possible scenario of homogeneously distributed perturbers of the same mass (Merritt 2013).Therefore, we test three typical values for the Coulomb logarithm, i.e. ln Λ = 5, 10, and 15.

Prior distributions of the input parameters
In total, a set of input parameters Θ consists of the DM scale radius R0,DM, the stellar scale radius R0, , the total DM mass MDM, and one or two parameters that characterise the assumed PBH mass function, i.e.Θ = {R0, , R0,DM, MDM, Mc} and Θ = {R0, , R0,DM, MDM, µ, σ} for the monochromatic and log-normal PBH mass function, respectively.Their prior distributions p(θi) are assumed to be independent of each other, so that the joint prior distribution p(Θ) in Eq. ( 1) is just the product of the individual priors: In order to reflect the increased variety of our observational sample (Table 1) compared to the one used by Zhu et al.
Thus, Eqs (2) and ( 13) can be used to determine the posterior distribution p(Θ | D) via Eq.( 1), that measures the likelihood of a simulated UFDG to meet the observations.

Monochromatic mass function
In this section, we describe the dynamical effects of monochromatic PBHs on the stellar population of the UFDGs.For this purpose, we perform four fiducial simulations varying the stellar scale radius R0, = {10, 20, 35, 50} pc, while fixing the mass of each PBH to Mc = 30 M , the total stellar mass to M = 10 3 M , the total DM mass to MDM = 10 9 M , and the DM scale radius to RDM,0 = 10 3 pc.We fix the Coulomb logarithm ln Λ = 15 for all of these simulations.
The resulting stellar mass density, ρ (r), and velocity dispersion, σ (r), profiles at the beginning and end of the simulation, as well as the temporal evolution of the deprojected half-mass radius R h, (t) are shown in Figure 2. The panels show that the central mass density of stars drops with time, while the central stellar velocity dispersion increases.This behaviour can be explained in terms of the initial temperature inversion of stars, i.e. the fact that their velocity dispersion exhibits a pronounced peak at some radius at the beginning of the integration (in our specific case, at ∼1 kpc).
In other words, one can think of the stars at the very centre as a dynamically "cool" sub-system (Binney 1980).As time progresses, gravitational encounters due to PBHs and other stars increase the central stellar velocity dispersion.As a result, the central stars get dynamically heated up and expand to larger radii.This phenomenon is well studied for the case of compact stellar systems lacking a central massive black hole, e.g.nuclear star clusters (Quinlan 1996).The expansion of the stellar population results in an enhancement of the stellar half-mass radius (right-hand panel in Figure 2).This phenomenon was analytically described by Brandt (2016): the stellar system slowly expands until it gets completely DM-dominated at all radii, while the half-mass radius grows with time as R h, ∝ t 0.5 .Our simulations yield a slightly less efficient heating rate for initial stellar scale radii R0, between 10 and 50 pc, as previously found by Zhu et al. (2018) for R0, up to 15 pc.The reason for this could be that two dynamical effects of the PBHs on the stellar population were implicitly neglected in the original derivation of the analytic expression by Brandt (2016): on the one hand, dynamical friction of the stars due to the PBHs cools the stellar population; on the other hand, the increase of the stellar velocity dispersion should gradually lower the heating efficiency by PBHs.If both effects are taken into account, heating would be, in principle, less efficient and the stellar population would expand at a lower pace than anticipated.In order to investigate whether these effects are significant or not, we calculate the complete heating-to-cooling ratio in the inner 10 2 pc at the initial time, t = 0, and final time, t = 12 Gyr.That is, we evaluate the fraction, for all of the four simulations as a function of radius, where D[(∆v ) 2 , D[(∆v ⊥ ) 2 ], and D[∆v ] (dynamical friction) are the diffusion coefficients for the stars that are assumed to follow a Maxwellian velocity distribution, X ≡ v/( √ 2σPBH) is the ratio between the typical stellar velocity v and the velocity dispersion σPBH of the PBHs, erf(X) is the error function, and the function G(X) is given by (Binney & Tremaine 2008): The velocity dispersion σPBH is a direct output of Phase-Flow as a function of radius, whereas for the former we assume that v = √ 3σ = √ 3σ (r), i.e. that it is given by their root-mean-square speed.We find that the heating-tocooling ratios are bigger than one in all simulations but they all drop on the course of time -sometimes drastically.In fact, the ratios drop from ∼14 (t = 0) to ∼11 (t = 12 Gyr) for R0, = 50 pc, from ∼100 to ∼11 for R0, = 35 pc, from ∼300 to ∼12 for R0, = 20 pc, and from ∼800 to ∼12 for R0, = 10 pc.In all cases, the ratios appear to be constant with radius.Remarkably, all of the final cooling rates amount to almost 10 per cent of the corresponding heating rates -even though they were initially significantly lower.Therefore, one has to take into account the cooling due to dynamical friction and the increase of the stellar velocity dispersion, in order to simulate the dynamical evolution reliably.Finally, the central stellar values (density and velocity dispersion) plotted in Figure 2 are fairly insensitive to the initial stellar scale radius R0, after 12 Gyr.The final half-mass radius range also shrinks, compared to the initial range, albeit less so than in the case of the central stellar quantities.In Figure 3 we show the corner plot of the posterior distributions for R0,DM, R0, , MDM, and Mc resulting from the MCMC analysis on the stellar velocity dispersion and halflight radius of Leo IV (cf.Zhu et al. 2018, figure 5).There, we use a stellar mass-to-light ratio of 2 M /L and ln Λ = 15.The posteriors for quantities which followed a log-uniform prior (R0,DM, R0, , and Mc) quickly develop a pronounced peak, whereas the posterior for MDM is dominated by our choice of a log-normal prior.The light and dark brown contours enclose 2σ and 1σ areas, respectively.Analogously, the filled areas below the curves indicate the 1σ interval with respect to their maxima.We emphasise that the corner plot we show for Leo IV is exemplary for the other UFDGs, as their posterior distributions have a similar shape.Most importantly for this work, all the UFDGs develop pronounced peaks in the posterior distributions for the PBH mass Mc.We take the respective 1σ interval around these peaks as the PBH mass region preferred by the respective UFDG.For Leo IV, a wide range, Mc ∈ [8, 162] M , is thus allowed.In Figure 4, we show the PBH mass regions preferred by each UFDG separately and also present the results for different values of the stellar mass-to-light ratio (1, 2, and 3 M /L ) and of the Coulomb logarithm (5, 10, and 15).The UFDGs are sorted from top to bottom from the lightest (faintest) to the heaviest (brightest) galaxy; from this, a clear tendency is inferrable: light (faint) UFDGs prefer light PBHs, whereas heavy (bright) UFDGs prefer heavy PBHs.This can be concluded for all stellar mass-to-light ratios and Coulomb logarithms.In fact, the intersection of all preferred mass intervals is always empty.Put differently, there exists no single value of the PBH mass that meets all preferences simultaneously, as there are always multiple UFDGs disfavouring it.This is a key result of this paper.In the figure, we also show the mass interval Mc ∈ [25, 100] M allowed by previous constraints (Carr et al. 2017).2For all choices of the stellar mass-to-light ratio and of the Coulomb logarithm, the UFDGs Draco II, Segue I, Tucana III, Triangulum II, Segue II, and Grus I are incompatible with the previous constraints.Considering the default values (stellar mass-tolight ratio = 2 M /L and ln Λ = 15), even the nine faintest UFDGs and Leo V turn out to be incompatible with that mass interval.When using the MCMC analysis in order to infer the preferred parameters for PBH-DM, we implicitly assume that the DM consists of PBHs, i.e. we do not question the very existence of PBHs but only the likelihood of the input parameters.In fact, given an UFDG, the inferred preferred ranges for its input parameters R0,DM, R0, , MDM, and Mc do not necessarily lead to the observed stellar half-mass radius R h, and velocity dispersion σ within the reported uncertainty (cf.Table 1).Instead, they could lead to values that are just not as bad as all the other values in the explored parameter space.Nevertheless, they might sit at the low-probability tails of Eq. ( 2).For instance, running a simulation (stellar mass-to-light ratio = 2 M /L , ln Λ = 15) of Leo IV with its likeliest input parameters, R0,DM = 103.16pc, R0, = 10 2.03 pc, MDM = 10 9.03 M , and Mc = 10 1.84 M , i.e. the maxima of the posterior distributions in Figure 3 (cf.Table 2), yields the stellar velocity dispersion σ ,S = 3.6 km s −1 and the deprojected half-mass radius R h, ,S = 181.4pc. 3 This result has to be compared with the observed values σ ,D = 3.3 ± 1.7 km s −1 and R h, ,D = (4/3)r h, ,D = 152 ± 16 pc (cf.Table 1).While the simulated and observed stellar velocity dispersions agree well with each other, the proposed half-mass radius is actually too large to meet the observed one within the given measurement uncertainty.We repeat this analysis for each individual UFDG and add the results to Table 2.For comparison, we also add the observed values adopted from Table 1.It can be seen that the proposed stellar velocity dispersions σ ,S of a large majority of UFDGs (25 out of 27) match the observed values within the respective measurement uncertainty, the exceptions being Reticulum II and Hercules, for which the simulations yield values that are slightly too large.In contrast, the proposed half-mass radii for 16 out of 27 UFDGs do not agree with the observed values within the measurement uncertainty.They are: Draco II, Segue II, Willman I, Reticulum II, Pisces II, Ursa Major II, Aquarius II, Coma Berenices, Carina II, Hydra II, Hydrus I, Leo IV, Canes Venatici II, Hercules, Boötes I, and Eridanus II.The simulations of all of them except Draco II output a value which is too large.This outcome constrains the monochromatic PBHs mass functions even more severely than the result from the MCMC analysis described in the previous paragraph.There, we concluded from the combined analysis of different UFDGs that there exists no single PBH mass that works for all of them at the same time, although we implicitly assumed that PBHs constitute DM in each UFDG.Here, in turn, we find that the majority of UFDGs conflicts this very assumption individually since PBHs cannot reproduce their stellar half-mass radii at all.Table 2. Preferred input parameters R 0,DM , R 0, , M DM , and Mc for each UFDG separately.They are evaluated as the 1σ intervals around the maxima of the respective posterior distributions.For each individual UFDG, a simulation has been performed with the given input parameters yielding a deprojected stellar half-mass radius and stellar velocity dispersion, reported as R h, ,S and σ ,S , respectively.For comparison, we also add the observed values, R h, ,D = (4/3)r h, ,D and σ ,D , adopted from Table 1.The PBH mass function is monochromatic and the stellar mass-to-light ratio and the Coulomb logarithm are 2 M /L and 15, respectively.

Log-normal mass function
Here, we investigate whether the constraints resulting from the MCMC with an extended, log-normal mass function for PBHs weaken those from the previous section.We restrict our analysis to a stellar mass-to-light ratio of 2 M /L and a Coulomb logarithm of 15.
In Figure 5, we present the resulting corner plot for Leo IV as a representative for the others.Note that we now deal with five input parameters, R0,DM, R0, , MDM, µ, and σ, instead of four since the log-normal mass function depends on two parameters.Notably, a comparison with the posterior distributions for R0,DM, R0, , and MDM that resulted from the monochromatic simulations (Figure 3, Table 2) shows no significant difference.In fact, the 1σ intervals around the respective maxima of the posterior distributions which are listed in Table 3 agree well with the previous one: log 10 (R0,DM/pc) = 3.14 +0.39 −0.30 (3.16 +0.39 −0.31 for the monochromatic case), log 10 (R0, /pc) = 2.03 +0.12 −0.26 (2.03 +0.11 −0.22 ), and log 10 (MDM/ M ) = 8.90 +0.58  −0.45 (9.03 +0.45 −0.59 ).Concerning the posterior distribution for µ, we have to bear in mind that it is, in general, not directly comparable to Mc, since they have different meanings.The latter is equivalent to the mean mass of the monochromatic mass distribution, whereas this is given by M = µ exp(−σ 2 /2) for the log-normal.At this point, we refer to Green (2016), who proposed to obtain a general constraint on an extended mass function by a replacement that, in our case (ln Λ = 15 and fPBH = 1), reads: For Leo IV, log 10 σ = −0.75+0.56 −0.17 is small, meaning that µ exp(+σ 2 /2) ≈ µ.Indeed, log 10 (µ/ M ) = 1.81 +0.39  −1.03 agrees well with log 10 (Mc/ M ) = 1.84 +0.37  −0.96 , confirming the proposed method by Green (2016).In Figure 6, we show the preferred input parameters for all UFDGs separately and add them to Table 3.We find no value of µ that is allowed by all UFDGs, analogously to what found in Section 3.1.Again, there is the tendency for light (faint) UFDGs to prefer light PBHs and heavy (bright) UFDGs to prefer heavy PBHs.The nine lightest UFDGs still constitute the critical sub-sample that by itself would conflict with the previous constraints (Carr et al. 2017).In turn, the constraints of Carr et al. (2017) agree well with the preferred ranges for σ of all our UFDGs.In fact, there is much scatter across the allowed range, suggesting that the constraint is fairly insensitive to the actual value of σ.If the result depends on µ exp(+σ 2 /2), as suggested by Eq. ( 23), this is reasonable as the contribution of σ would be exponentially suppressed for all UFGDs.Finally, we also investigate whether the likeliest values for the input parameters can be used to reconstruct the observed properties of the stellar populations.For this purpose, we proceed in the same way as in Section 3.1 and add the resulting stellar half-mass radii R h, ,S and stellar velocity dispersions σ ,S to Table 3. Again, the observed stellar velocity dispersions σ ,D can be well reproduced.Indeed, all UFDGs agree within the given measurement uncertainty.In turn, the half-mass radii of 16 UFDGs are incompatible with the observed, R h, ,D .From the analysis above, we conclude that the log-normal mass function does not weaken the constraints we found for the monochromatic mass function.In fact, the width of the distribution seems to have little or no effect on the dynamical constraints.5).The PBH mass function is log-normal and the stellar mass-to-light ratio and the Coulomb logarithm are 2 M /L and 15, respectively.Note that the prior distributions for R 0,DM , R 0, , µ, and σ are log-uniform, whereas that for M DM is log-normal.The light and dark orange contours enclose the 2σ and 1σ areas, respectively.Analogously, the filled areas below the curves indicate the 1σ interval with respect to their maxima.Table 3.The likeliest input parameters R 0,DM , R 0, , M DM , µ, and σ from the MCMC with a log-normal mass function (cf.Table 2 for the monochromatic case).They are evaluated as the 1σ intervals around the maxima of the respective posterior distributions.Again, we performed for each individual UFDG a simulation with these values resulting in the listed values for R h, ,S and σ ,S .For comparison, we also add the observed values, R h, ,D = (4/3)r h, ,D and σ ,D , adopted from

DISCUSSION AND CONCLUSIONS
In this work, we investigated the dynamical effect of O(1-100) M PBH-DM on the stellar population of UFDGs by means of Fokker-Planck simulations.We showed that PBHs heat up their stellar system by means of two-body scatterings that result in an increase of the stellar velocity dispersion and in the growth of the stellar half-mass radius.We inferred the posterior probability distributions for certain input parameters to reconstruct the observed stellar velocity dispersions and half-mass radii of UFDGs by means of 2×10 5 MCMC runs per galaxy.In particular, we investigated the effect of different PBH mass functions (monochromatic and log-normal) on the dynamical evolution of stars inhabiting DM-dominated UFDGs.
As a result, we are able to make three statements each of which independently challenges the PBH-DM scenario: (i) Concerning the monochromatic mass function, we found that the lowest-mass (faintest) UFDGs constitute a critical sub-sample that prefers a PBH mass range Mc ∼ O(1) M ; this range is compatible with some previous dynamical constraints (Brandt 2016;Zhu et al. 2018) but incompatible with other constraints (Carr et al. 2017).Assuming a log-normal PBH mass function does not resolve this conflict.(ii) Even disregarding the previous point, the combined analysis of all 27 UFDGs in our sample leaves no open mass window for Mc (the mass for the monochromatic mass function) and µ (the scale mass for the log-normal mass function) in which PBH-DM could simultaneously produce the observed stellar properties for all galaxies.(iii) In a majority of UFDGs, PBH-DM would puff up the stellar systems to half-mass radii that are too large to meet the observed values (even accounting for the uncertainty in the observables).
In the following, we discuss the possible shortcomings of our work.First of all, the stellar observables of several among the lightest UFDGs suffer from poor statistics (as mentioned in Table 1) due to the small number of detectable stars.This could in principle weaken our point (i) above; however, we stress that it would be extremely unlikely that all the stellar observables in this sub-sample would shift our results in the same direction.Similarly, we note that the very nature of some of the lightest UFDGs in our sample is less clear, e.g. that of Tucana III (Simon et al. 2017) and Triangulum II (Kirby et al. 2017).To the present, it cannot be ruled out that those are actually baryon-dominated star clusters instead of DM-dominated UFDGs.If future observations corroborate this scenario, the respective systems must not be included in our sample in order to draw meaningful conclusions on PBHs as DM.Our results can be easily updated to such situation by simply removing the respective rows in Tables 2 and 3 and Figures 4 and 6.However, we highlight that none of our three main findings would change unless a majority of UFDGs in the critical sub-sample turns out to be spurious.Furthermore, we emphasise that we only studied the scenario in which PBHs constitute all of the DM ("PBH-DM").
For models in which DM consists only partially of PBHs, the three statements might no longer be valid.In particular, our method is inadequate to infer any constraint on the possibility that DM is complementarily composed of PBHs and a collisionless component, as our methodology does not allow us to implement the latter.
We additionally neglected the possible influence of gas and a central, intermediate-mass black hole on the dynamics of UFDGs.The former assumption is well-justified as the observed UFDGs show no significant gas content [cf.Simon (2019)].On the other hand, a massive black hole could inhabit the core of at least some among our galaxies (see, e.g. the recent observation by Baldassare et al. 2017 of a more massive dwarf).However, the mass of a central black hole in such low-mass systems is expected to be very small ( 10 4 M ) if we extrapolate the commonly accepted scaling relations (Shankar et al. 2016) to the lowest mass end.Such black holes are expected to have a non-negligible dynamical influence only at very small scales ( 1 pc), where a small cusp could develop (Bahcall & Wolf 1976).On the other hand, such black holes may widely oscillate around the centre of the system as they are only 10 4 times heavier than the field population.Such wandering could in principle dynamically heat and puff up the system even more.
In our analysis, we further assumed that the UFDGs are spherically symmetric and in dynamical equilibrium.In fact, only a handful of them exhibit significant flattening and/or signs of rotation (McConnachie 2012).Tidal stirring and mass stripping from the primary halo could at first make both the stellar and dark matter distribution more elongated; however, dynamical relaxation is well known to make the system round over a sufficiently long time-scale even if the system was initially maximally flat, as was repeatedly shown by numerical simulations of both galaxy clusters (Mastropietro et al. 2005) and Milky Way-like halos (Tomozeiu et al. 2016;Buck et al. 2018).The assumption of spherical symmetry made throughout this work thus appears reasonable.Concerning the dynamical state of the UFDGs, it is expected that phenomena which would put the UFDGs out of equilibrium, e.g.tidal stirring, tendentially lead to extra dynamical heating of the stars (Mayer 2010;Tomozeiu et al. 2016).Therefore, the stellar systems would puff up further and exhibit increased velocity dispersion and halfmass radii, thus making our findings on the evolution of both quantities conservative.
Having discussed some possible shortcomings of our method, we highlight that our results join the class of dynamical constraints that severely challenge the PBH-DM scenario.The PBH mass window previously left open was already in conflict with other constraints (Carr et al. 2017).Here, we closed this window by improving previously applied methods in several ways.In particular, we were also able to show that the increase of the stellar velocity and the dynamical friction due to PBHs actually lower the heating efficiency in a non-negligible way, making our results more precise than previous ones.
In light of these constraints, the PBH-DM scenario appears unlikely to be realised in our Universe, and the very nature of DM remains elusive.

Figure 1 .
Figure 1.Binning scheme for a log-normal PBH mass function with µ = 20 M and σ = 0.7.The total DM mass is set to M DM = 10 9 M .The 1 per cent of the mass in the tails of the distribution is neglected.Thus, the N = 50 bins are limited by the thresholds M th,0 ≈ 3.3 M and M th,1 ≈ 121.4 M and in total make up 99 per cent of the total DM mass.Each bin corresponds to one PBH-DM component in the simulation with a PBH mass approximated with the average PBH mass within the bin, cf.Eq. (11).

Figure 2 .
Figure 2. Stellar mass density profile (left-hand panel), stellar velocity dispersion profile (central panel), and temporal evolution of the deprojected half-mass radius R h, (right-hand panel) for four runs with a monochromatic PBH mass function (Mc = 30 M ), a total DM mass of M DM = 10 9 M , a total stellar mass of M = 10 3 M , and a DM scale radius R DM,0 = 10 3 pc, which differ only in the initial stellar scale radius: R 0, = 10 (green), 20 (red), 35 (blue), and 50 pc (yellow).The left-hand and central panels show the stellar profiles at the beginning (t = 0; dashed lines) and at the end of the simulation (t = 12 Gyr; solid lines).The grey dotted and dashed lines in the right-hand panel show the reference growth rates R h, ∝ t 0.4 and R h, ∝ t 0.5 , respectively.

Figure 3 .Figure 4 .
Figure 3. Corner plot of the posterior distributions for R 0,DM , R 0, , M DM , and Mc, based on the stellar velocity dispersion and halflight radius of Leo IV (cf.Zhu et al. 2018, figure 5).The PBH mass function is monochromatic and the stellar mass-to-light ratio and the Coulomb logarithm are 2 M /L and 15, respectively.Note that the prior distributions for R 0,DM , R 0, , and Mc are log-uniform, whereas that for M DM is log-normal.The light and dark brown contours enclose the 2σ and 1σ areas, respectively.Analogously, the filled areas below the curves indicate the 1σ interval with respect to their maxima.

Figure 5 .
Figure 5. Corner plot of the posterior distributions for R 0,DM , R 0, , M DM , µ, and σ based on the stellar velocity dispersion and half-mass radius of Leo IV (cf.Zhu et al. 2018, figure5).The PBH mass function is log-normal and the stellar mass-to-light ratio and the Coulomb logarithm are 2 M /L and 15, respectively.Note that the prior distributions for R 0,DM , R 0, , µ, and σ are log-uniform, whereas that for M DM is log-normal.The light and dark orange contours enclose the 2σ and 1σ areas, respectively.Analogously, the filled areas below the curves indicate the 1σ interval with respect to their maxima.

Figure 6 .
Figure 6.Preferred 1σ intervals of the input parameters for the log-normal PBH mass function.The stellar mass-to-light ratio and the Coulomb logarithm are 2 M /L and 15, respectively.The fourth column (log 10 µ) has to be compared with the corresponding column (2 M /L , ln Λ = 15) of Figure 4. Still, there exists no single value of the PBH mass that meets all preferences simultaneously since the stellar kinematics appear to be fairly insensitive to the width σ of the mass function.The grey-shaded areas indicate the µ-window [25, 100] M and σ-window [0.0, 1.0] allowed by previous constraints (Carr et al. 2017).