© The Authors 2024

1. Introduction

The black hole (BH) mass spectrum is commonly divided into three categories, depending on the BH mass (M_BH), namely stellar mass black holes (sMBHs) with M_BH ≲ 100 M_⊙, supermassive black holes (SMBHs) with M_BH ≳ 10⁵ M_⊙, and intermediate-mass black holes (IMBHs), comprising everything in between (10² M_⊙ ≲ M_BH ≲ 10⁵ M_⊙; Greene et al. 2020). However, while sMBHs and SMBHs have been widely detected, the IMBH regime remains unpopulated, with very few detections based on gravitational waves (Abbott et al. 2020) and several controversial cases (Gebhardt et al. 2005; van der Marel & Anderson 2010). This observational lack is somewhat remarkable considering that the Universe must have been, at some point, inhabited by IMBHs, as the observed SMBHs are believed to have grown from less massive ones. Additionally, the extremely massive SMBHs that have been discovered in young quasars at high redshifts (Bañados et al. 2018) are even more challenging to explain in the absence of massive seeds since it is unclear how SMBHs could have grown so fast in time otherwise.

Based on extrapolations from stellar mass–BH mass relations, the prevailing notion is that dwarf galaxies and globular clusters (GCs) are prime locations to search for IMBHs. However, despite considerable efforts dedicated to this quest, detections remain elusive (Strader et al. 2012; Lützgendorf et al. 2016; Della Croce et al. 2024). In these low-mass systems, significant challenges lie, for instance, in requesting high spatial resolution to collect data from a large number of kinematic tracers close to the BH, in the accurate identification of the system’s center (Maccarone et al. 2005; Giersz et al. 2015; de Vita et al. 2017; Pechetti et al. 2024), or in the assumption of correct stellar mass-to-light ratios. Nonetheless, despite the observational and theoretical obstacles, there has been a recent claim regarding the discovery of an exceptionally massive BH at the center of the local dwarf spheroidal galaxy (dSph) Leo I (Bustamante-Rosell et al. 2021, hereafter BR21). According to BR21, who employed the orbit-based Schwarzschild (1979) method and a frequentist model–data comparison, the 1σ BH mass is M_BH = (3.3 ± 2)×10⁶ M_⊙, comparable to that of Sagittarius A^*. This finding is exceptional for several reasons: (i) it is an actual detection, not just an upper limit, with the no-BH case excluded at 95% significance; (ii) the mass of the BH is comparable with the stellar mass of the entire galaxy; and (iii) the BH would fall into the category of SMBHs and not IMBHs.

This recent discovery poses a challenge, as the BH in Leo I deviates from any scaling relation by at least two orders of magnitude (Greene et al. 2020). This exceptional deviation calls for an explanation, and therefore further investigations into Leo I to gain a deeper understanding of its unique properties are warranted. For instance, Pacucci et al. (2023) recently suggested that Leo I might have been more massive in the past, with pericentric passages close to the Milky Way (MW) potentially stripping up to 90% of its original mass. However, while their simulations reproduce Leo I’s current position and velocity dispersion, they require very strict and improbable conditions: two pericentric passages – a scenario that is unlikely for the galaxy (Sohn et al. 2013) – that are systematically close to the lower bound of the pericenter as estimated from Gaia Early Data Release 3 data (Pace et al. 2022). In this Letter we tackle the problem in a different manner, and we argue that the structure and kinematics of Leo I are, instead, perfectly consistent with a no-BH scenario and that the BH, if present, is at most of intermediate mass. In our work, we fit the same dataset as in BR21 but using different dynamical models, based on distribution functions (DFs) that depend on actions, in a fully Bayesian framework.

Leo I appears as a very regular, flattened structure on the plane of the sky, with an axis ratio of 0.69 ± 0.01, similar to that of other classical dSphs (Muñoz et al. 2018). It is among the brightest and most massive dSphs in the MW, with a stellar mass of 5.5 × 10⁶ M_⊙ (McConnachie 2012). Its heliocentric distance is D = 256.7 ± 13.3 kpc (Méndez et al. 2002; Bellazzini et al. 2004; Held et al. 2010; Pacucci et al. 2023), which designates the galaxy as the most remote dSph within the MW. Throughout this work we adopt D = 256.7 kpc.

This Letter is structured as follows: In Sects. 2 and 3 we introduce the used dataset and models, respectively. Section 4 describes the fitting procedure, while in Sect. 5 we present our results and draw our conclusions.

2. Data

The dataset used in this work consists of the system’s surface brightness profile, the central line-of-sight velocity distributions (LOSVDs), and the galaxy outer velocity dispersion profile.

Photometric sample. We adopted the surface brightness profile presented in BR21, which was derived using publicly available Sloan Digital Sky Survey (SDSS) g-band imaging (Data Release 12; Alam et al. 2015). This profile extends from 4 arcsec – where, as pointed out by BR21, it shows a mild cusp – to 20 arcmin (≃1.5 kpc at the adopted distance). The profile, shown in the right-hand panel of Fig. 1, thus consists of N_μ = 26 radial bins, the j-th of which has an average distance R_j, associated V-band magnitude μ_V, j, and error Δμ_V, j.

Fig. 1. Comparison between data and models. Bottom left: kinematic sample from M08. In blue (gold), we show stars included in (excluded from) our analysis. Top left: derived los velocity dispersion profile. The bins’ edges are indicated by vertical dashed lines in the bottom-left panel. The red points represent the los velocity dispersions computed from the LOSVDs shown in Fig. 2. Middle panel: marginalized velocity distribution. Top right: surface brightness profile from BR21. Bottom right: model residuals compared to the data. The light (dark) blue and red bands represent the models’ 1σ (3σ) confidence intervals, while the solid curves are median models.

Inner kinematic sample. The kinematic sample used in the central galaxy encompasses the galaxy’s LOSVDs from BR21¹. These distributions were computed from observations collected with the VIRUS-W spectrograph of two overlapping 1.75  ×  0.92 arcmin² rectangular patches in the galaxy’s central regions. Due to spectrograph fiber size limitations, individual stellar spectrum extraction was unfeasible, so integrated-light kinematics was employed. To mitigate crowding effects, the fibers were grouped into N_R = 23 sectors, each containing from 5 to 41 fibres, distributed across a polar grid centered on the galaxy center (see BR21 for details). Figure 2 shows the spatial distribution of the sectors (top right) along with each LOSVD. Here, the i-th LOSVD has a number N_v, i of velocity bins, and the k-th bin has a probability ℒ_k, i and error Δℒ_k, i.

Fig. 2. LOSVDs from BR21 (colored histograms) compared to those from the median models (red curves). The columns, from left to right, show LOSVDs computed at progressively larger radii. Rows, from top to bottom, show LOSVDs computed at increasing position angles. The bottom row (in orange) displays LOSVDs from BR21 weighted along the position angle of each radius. The top-right part of the figure shows the spatial distribution of the LOSVDs, color-coded for the offset velocity.

Outer kinematic sample. For the galaxy’s outer parts, we employed a line-of-sight (los) velocity dispersion profile based on the kinematic sample from Mateo et al. (2008, hereafter M08). The profile consists of N_σlos = 8 radial bins, the k-th of which has an average distance R_k, velocity dispersion σ_los, k, and error Δσ_los, k. Figure 1 shows the sample from M08 and the derived los velocity dispersion profile. Details on how the profile was computed are given in Appendix B.

3. Models

Leo I is represented as a multi-component galaxy with stars, dark matter (DM), and a possible BH in its center. The phase-space distribution of stars and DM is described with a DF that depends on the action integrals, J,

(1)

where

(2)

with the DM labelled as i = dm and the stars as i = ⋆ (Posti et al. 2015; Cole & Binney 2017). DFs that depend on actions have a history of successful application: Binney & Vasiliev (2023, 2024) applied chemo-dynamical models based on f(J) DFs to model the phase-space distribution of the MW; Pascale et al. (2018) first applied these DFs to the Fornax dSph, showing that its stellar kinematics is well described only by a cored DM density distribution (see also Pascale et al. 2019); and Della Croce et al. (2024) employed f(J) models to place a very tight upper limit on the mass of a putative IMBH at the center of the GC 47 Tucanae.

DF (1) provides the flexibility to characterize components whose density distributions are well described by a double power-law model (with the optional presence of a central density core) and are truncated in the outer parts. The velocity distribution varies with the distance from the center, shifting from one velocity bias (tangential, isotropic, or radial) in the center to another in the outer parts. The DF (Eq. (1)) is normalized to the target component total mass, M_i, via the normalization constant f_i; J_c, i sets the size of the central component’s cored region in action space and marks the transition between the cored region and the inner power law, with slope Γ_i. We note that when J_c, i = 0, the model has no density core. Around J_i, the outer power law of index B_i − Γ_i sets in. This transition happens with a sharpness described by η_i > 0. The DF is quickly truncated around ≃J_t, i, with a strength given by α_i ≥ 0. The parameters 0 < h_r, i < 3, 0 < h_z, i < 3 attribute different weights to radial and vertical orbits, thus regulating the inner velocity distribution of the models. The same task is achieved with 0 < g_r, i < 3, 0 < g_z, i < 3, but in the model’s outer parts. The inner and outer slopes and the truncation of the DF in the action space lead to similar behaviors in the physical space.

We accounted for the self-gravity of both stars and DM by solving the Poisson equation

(3)

where i ∈ { ⋆ ,dm}, x is the position vector and v is the velocity. The total gravitational potential is

(4)

with

(5)

describing a third, point-like, fixed, contribution added to capture the optional presence of a central BH (G is the gravitational constant and vector and r ≡ ||x||).

In the formalism of Eq. (2), the models can be flattened along the z-axis by adjusting the weights appearing in the g_i(J) and h_i(J) functions. Here, however, we limit ourselves to spherical models requiring

(6)

and thus

(7)

In other words, the DFs depend on actions only through the angular momentum magnitude L = |J_ϕ|+J_z. Without loss of generality, in the following we assume that z is the los and that R is the distance, on the plane of the sky, from the galaxy center.

Once the DF is specified, appropriate integrations yield structural or kinematic quantities that can be compared with data or used to make predictions. Key examples include the stellar projected density,

(8)

the stellar LOSVD,

(9)

and its second moment (i.e., the los velocity dispersion),

(10)

where v_⊥ is the component on the plane of the sky of the velocity vector, v, and is perpendicular to the los velocity, v_los (where the average los velocity is null since the models are nonrotating). We relied on the AGAMA² (Vasiliev 2019) software library to solve all these integrals and to compute the models’ total potential.

4. Model fitting

The log-likelihood of the model, lnℒ, given the available set of data, 𝒟, is

(11)

Explicitly,

(12)

(13)

and

(14)

Each term addresses a distinct dataset (see Sect. 2): Eq. (12) models the galaxy’s surface brightness profile using the DF-based model (Eq. (8)), suitably scaled to the V-band surface brightness via the nuisance parameter, c (which is proportional to the V-band stellar mass-to-light ratio); Eq. (13) fits the los velocity dispersion profile with the σ_los derived from Eq. (10); Eq. (14) the fit to the central LOSVDs.

A single galaxy model is entirely determined by 17 free parameters:

• M_⋆, J_⋆, J_c, ⋆, h_z, ⋆, g_z, ⋆, Γ_⋆, B_⋆, and η_⋆, which belong to the stellar DF (1), with i = ⋆.

• ], J_dm, J_c, dm, g_z, dm, η_dm, Γ_dm, and B_dm of the DM DF (1), with i=dm.

• The BH mass, M_BH.

• The nuisance parameter, c.

To decrease the number of free parameters, we fixed

• α_⋆ = 0 in the stellar DF since data do not require an exponential cutoff.

• α_dm = 2 and J_t,dm = 200 kpc km s⁻¹ in the DM DF. Here, since the outer halo density cannot be constrained by the available dataset, it is meaningless to fit the halo truncation radius. By truncating the halo, we also ensure that the total DM mass is a well-defined quantity. For spherical models, the exact values of J_t, dm and α_dm do not matter as long as the halo is truncated beyond the radial extent of the stellar component.

• h_z, dm = g_z, dm since the inner and outer anisotropy of the DM halo is an unconstrained quantity in the present dataset.

We ran a Markov Chain Monte Carlo (MCMC) procedure to sample from the posterior distribution. Further details are given in Appendix C, and the prior and the 1σ and 3σ confidence intervals on the models’ parameters resulting from the analysis are listed in Table A.1.

5. Results and conclusions

Here we present and discuss the key results of our study. Figure 3 shows the one-dimensional, marginalized, posterior distribution on log M_BH. Over the explored prior (see Table A.1), the posterior is almost uniform toward low masses and slightly increases around log M_BH ≃ 5.5, whereas it has a sharp truncation around ∼10⁶ M_⊙. This strongly suggests that there is insufficient statistical evidence to assert the detection of a BH. Our inference can only go so far in establishing an upper limit on the BH mass: 6.76 × 10⁵ M_⊙ at 3σ. Our 3σ upper limit is lower than the 1σ lower limit reported by BR21 (≃1.3 × 10⁶ M_⊙) but, as can be seen in Fig. 3, is marginally consistent with the 3σ lower bound derived from the reference model of BR21 (see the vertical band in our Fig. 3). We note that, in addition to their reference model, BR21 present other three models, and for all three they infer a 1σ lower limit on M_BH that is higher than 10⁶ M_⊙. However, upon careful examination of these models, it turns out that they are consistent with a no-BH scenario at the 3σ level. We stress that BR21 do not explicitly report the 3σ intervals on their models’ free parameters, so we computed them directly from their Fig. 11, and in particular from the panels that plot χ² as a function of the BH mass. Consistent with the frequentist approach used in BR21, we computed the 3σ interval on the BH mass considering the χ² of all the models with Δχ² = 3, where Δχ² is the χ² difference from the model with the minimum χ². The model with the minimum χ² corresponds to the best-fit model.

Fig. 3. Posterior distribution on the logarithm of the BH mass, log MBH (represented by the blue curve). The blue horizontal bar shows the 3σ confidence interval. The posterior is compared with the best model and with the 1σ and 3σ estimates from BR21 (shown by the solid orange line and bands).

When compared to the spatial extent of the LOSVD dataset, the inferred BH mass translates to an upper limit on the radius of influence R_infl ≡ GM_BH/σ_los (σ_los is los velocity dispersion) of the BH that is comparable with the average distance of the innermost kinematic sectors (≃15.35 arcsec = 19 pc). In the context of BR21’s study, the inferred BH mass instead implies a well-defined R_infl that is comparable to the distance of the outer LOSVD sectors. Therefore, in our case, the kinematics of the inner regions effectively rules out the existence of a SMBH rather than supporting its presence.

In Figs. 1 and 2, we provide an illustration of the agreement between the data and models. Figure 1 shows the models’ predictions for the los velocity dispersion and the surface brightness profiles compared to observations. Figure 2 depicts the LOSVDs from BR21 and the analogs from the models. In both cases, the models accurately fit the data.

Based on scaling relations (e.g., Fig. 3 from Greene et al. 2020), for a galaxy similar to Leo I, with a velocity dispersion between 10 and 12 km s⁻¹, one would expect an IMBH of, at most, 10⁴ M_⊙, or equivalently μ ≡ M_BH/M_⋆ = 10⁻⁴ − 10⁻³. Our upper limits on the BH mass and μ are in good agreement with these values. In terms of μ, we measure a median μ = 0.0013, with a 3σ upper limit at 0.033. As previously explained, these values represent the lower bounds that can be imposed by the kinematic dataset of LOSVDs.

The use of integral field spectroscopy to look for IMBHs in dense stellar systems has been extensively debated in the past. In GCs, it is believed, for instance, that the methodology may introduce biases, as the collected spectra can be dominated by a few bright stars rather than sampling the underlying stellar distribution. NGC 6388 is emblematic in this respect: the spectra from individual stars sample a central velocity dispersion ≃10 km s⁻¹, while values from integral field spectroscopy indicate a dispersion as high as 25 km s⁻¹ (Lützgendorf et al. 2011), which ends up being interpreted as a signature of massive BHs. The case of Leo I is undoubtedly more enigmatic, and it is highly challenging to conclusively pinpoint the reason behind the discrepancy between our results and those from BR21. In contrast to the case of NGC 6388, where two datasets yield different results, we employed the same dataset of LOSVDs as BR21, but with different fitting algorithm and models. Thus, the use of integral field spectroscopy cannot be the reason for these differences.

We believe that the differences in the results stem from an interplay of factors. As pointed out by BR21 (see also our Fig. 1), the central los velocity dispersion of Leo I inferred from the LOSVDs rises, showing a mild cusp with respect to the central los velocity that can be measured from individual star velocities, aligning with the findings of Lanzoni et al. (2013). While the application of integral field spectroscopy introduces a potential bias toward favoring the presence of a massive BH, the integration of more generalized models allows for improved marginalization over key stellar properties, such as stellar anisotropy or stellar and DM distributions. These models effectively mitigate the initial bias, thereby contributing to a more nuanced interpretation of the findings. A potential resolution of this issue may be to use individual star velocities that sample the region given by the influence radius of the BH. Nevertheless, the considerable distance of Leo I poses a significant challenge. Assuming the LOSVD shape is crucial for detecting a BH, Amorisco & Evans (2012) estimate that a sample of at least 100 los velocities with an error smaller than 0.2σ_los is required to recover any velocity distribution. In the case of Leo I, this implies that more than 100 stars with a 2 km s⁻¹ error in the los velocity, confined to a region smaller than 15 arcseconds, are needed to infer a BH with a mass of a few 10⁵ M_⊙, as constrained by our upper limit.

Acknowledgments

We thank the anonymous referee for the useful comments and suggestions, which helped improve the quality of the work. We are very thankful to Maria Jose Bustamante-Rosell for sharing data. This paper is supported by the Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC), project funded by European Union – NextGenerationEU – and National Recovery and Resilience Plan (NRRP) – Mission 4 Component 2 within the activities of Spoke 3 (Astrophysics and Cosmos Observations). The research activities described in this paper have been co-funded by the European Union – NextGeneration EU within PRIN 2022 project n.20229YBSAN – Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch.

References

Appendix A: Additional table

Appendix B: Derivation of the los velocity dispersion profile

The sample of radial velocities from M08 used to compute the galaxy’s los velocity dispersion profile consists of 387 stars in the region of the sky occupied by Leo I. They were obtained with the Hectochelle multi-object echelle spectrograph on the MMT. To keep our dataset as similar as possible to that of BR21, we removed all stars inside a 90.3 arcsec radius, which corresponds to the radial extent of the outer LOSVD sectors. Also, we removed all stars flagged as possible binaries or false positives by M08 (i.e., observations in which the velocity of the sky could have been measured instead of the target star velocity; see M08 for details). We performed an iterative 3σ clip on the remaining stars to remove further foreground stars. As discussed by M08, the distinction between members and foreground stars is very clear, that is to say, the systemic velocity of Leo I is very different from the average velocity of foreground stars in that direction (Δv ≃ 250 km s⁻¹). The final sample consists of 288 stars with a median error of 2 km s⁻¹.

Using these stars, we constructed the galaxy’s los velocity dispersion profile. Each radial bin contains 35 stars except the last bin, which contains 43 stars. We adopted the method from Pryor & Meylan (1993) to compute the galaxy velocity dispersion. The likelihood of having a velocity dispersion σ_los, k in the k-th radial bin is

(B.1)

In the above equation, v_k, i and δv_k, i are the i-th radial velocity and corresponding error, σ_los, k and v_los, k are the velocity dispersion and average velocity of the bin (i.e., the models’ free parameters), and the sum extends over the stars in the bin (N_σlos, k = 35 or 43). We ran an MCMC procedure very similar to that of Appendix C to sample from the posterior. We took as a measure of the velocity dispersion of the k-th bin the median value of the marginalized, one-dimensional, posterior distribution and, as the error, the average between the 84th-50th percentile difference and 50th-16th percentile difference.

Appendix C: MCMC fitting procedure

In the MCMC fit of Section 4, to sample from the posterior, we used 150 chains, each evolved for approximately 20000 steps. We combined the differential evolution proposal from Nelson et al. (2014) and the snooker proposal from ter Braak & Vrugt (2008), as implemented in the software library EMCEE (Foreman-Mackey et al. 2013). During each MCMC iteration, we randomly selected one of the two proposals, assigning a 60% probability to the differential evolution proposal with a stretch factor of , where n = 17 (the dimension of the parameter space), and a 40% probability to the snooker proposal with a 0.85 stretch factor. Empirically, this combination improves MCMC performance by reducing chain autocorrelation and increasing acceptance rates.

Figure C.1 shows the one- and two-dimensional posterior distributions over the models’ free parameters. A zoomed-in view of the posterior distribution over log M_BH is shown in Fig. 3. All walkers except one, which did not reach convergence, were used to build the posterior and evaluate all confidence intervals. The initial 5000 steps from each walker were discarded as burn-in, ensuring a low autocorrelation and model convergence. We further applied a thinning of 50 steps, approximately on the order of the chains’ autocorrelations length. The 1σ confidence intervals over the model parameters or any derived quantity are estimated as the 16th, 50th, and 84th percentiles of the corresponding distributions, and the 3σ confidence intervals as the 0.15th, 50th, and 99.85th percentiles.

Fig. C.1. Two- and one-dimensional marginalized posterior distributions over the model free parameters (see Table A.1). The dashed orange curves in the one-dimensional posterior distributions mark the 16th, 50th, and 84th percentiles, which were used to evaluate the 1σ confidence intervals of the corresponding parameters.

All Tables

All Figures

Fig. 1. Comparison between data and models. Bottom left: kinematic sample from M08. In blue (gold), we show stars included in (excluded from) our analysis. Top left: derived los velocity dispersion profile. The bins’ edges are indicated by vertical dashed lines in the bottom-left panel. The red points represent the los velocity dispersions computed from the LOSVDs shown in Fig. 2. Middle panel: marginalized velocity distribution. Top right: surface brightness profile from BR21. Bottom right: model residuals compared to the data. The light (dark) blue and red bands represent the models’ 1σ (3σ) confidence intervals, while the solid curves are median models.

In the text

Fig. 2. LOSVDs from BR21 (colored histograms) compared to those from the median models (red curves). The columns, from left to right, show LOSVDs computed at progressively larger radii. Rows, from top to bottom, show LOSVDs computed at increasing position angles. The bottom row (in orange) displays LOSVDs from BR21 weighted along the position angle of each radius. The top-right part of the figure shows the spatial distribution of the LOSVDs, color-coded for the offset velocity.

In the text

	Fig. 3. Posterior distribution on the logarithm of the BH mass, log MBH (represented by the blue curve). The blue horizontal bar shows the 3σ confidence interval. The posterior is compared with the best model and with the 1σ and 3σ estimates from BR21 (shown by the solid orange line and bands).
In the text

	Fig. C.1. Two- and one-dimensional marginalized posterior distributions over the model free parameters (see Table A.1). The dashed orange curves in the one-dimensional posterior distributions mark the 16th, 50th, and 84th percentiles, which were used to evaluate the 1σ confidence intervals of the corresponding parameters.
In the text