# Detection of a 100,000 M black hole in M31's Most Massive Globular Cluster: A Tidally Stripped Nucleus

Renuka Pechetti, Anil Seth, Sebastian Kamann, Nelson Caldwell, Jay Strader, Mark den Brok, Nora Luetzgendorf, Nadine Neumayer, Karina Voggel

► To cite this version:

Renuka Pechetti, Anil Seth, Sebastian Kamann, Nelson Caldwell, Jay Strader, et al.. Detection of a 100,000 M black hole in M31's Most Massive Globular Cluster: A Tidally Stripped Nucleus. The Astrophysical Journal, 2022, 924 (2), pp.48. 10.3847/1538-4357/ac339f. hal-03542062

HAL Id: hal-03542062

https://hal.science/hal-03542062

Submitted on 25 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Detection of a 100,000 $M_{\odot}$ black hole in M31's Most Massive Globular Cluster: A Tidally Stripped Nucleus

Renuka Pechetti<fnref n="1" />, Anil Seth<fnref n="2" />, Sebastian Kamann<fnref n="1" />, Nelson Caldwell<fnref n="3" />, Jay Strader⁴, Mark den Brok<fnref n="5" />, Nora Luetzgendorf<fnref n="6" />, Nadine Neumayer<fnref n="7" />, and Karina Voggel<fnref n="8" />

Received 2021 August 18; revised 2021 October 25; accepted 2021 October 25; published 2022 January 11

## Abstract

We investigate the presence of a central black hole (BH) in B023-G078, M31's most massive globular cluster. We present high-resolution, adaptive-optics assisted, integral-field spectroscopic kinematics from Gemini/NIFS that show a strong rotation (~20 km s⁻<fnref n="1" />) and a velocity dispersion rise toward the center (37 km s⁻<fnref n="1" />). We combine the kinematic data with a mass model based on a two-component fit to HST ACS/HRC data of the cluster to estimate the mass of a putative BH. Our dynamical modeling suggests a >3σ detection of a BH component of 9.1<sup>+2.6</sup><sub>-2.8</sub> × 10⁴ $M_{\odot}$ (1σ uncertainties). The inferred stellar mass of the cluster is 6.22<sup>+0.03</sup><sub>-0.05</sub> × 10⁶ $M_{\odot}$, consistent with previous estimates, thus the BH makes up 1.5% of its mass. We examine whether the observed kinematics are caused by a collection of stellar mass BHs by modeling an extended dark mass as a Plummer profile. The upper limit on the size scale of the extended mass is 0.56 pc (95% confidence), which does not rule out an extended mass. There is compelling evidence that B023-G078 is the tidally stripped nucleus of a galaxy with a stellar mass >10⁹ $M_{\odot}$, including its high-mass, two-component luminosity profile, color, metallicity gradient, and spread in metallicity. Given the emerging evidence that the central BH occupation fraction of >10⁹ $M_{\odot}$ galaxies is high, the most plausible interpretation of the kinematic data is that B023-G078 hosts a central BH. This makes it the strongest BH detection in a lower-mass (<10⁷ $M_{\odot}$) stripped nucleus, and one of the few dynamically detected intermediate-mass BHs.

*Unified Astronomy Thesaurus concepts: Andromeda Galaxy (39); Globular star clusters (656); Intermediate-mass black holes (816); Stellar kinematics (1608); Local Group (929)*

## 1. Introduction

Intermediate-mass black holes (IMBHs) are hypothesized to exist in the mass range between stellar-mass black holes (≈100 $M_{\odot}$) and super-massive black holes (SMBHs; ≳10⁵ $M_{\odot}$). Some models of SMBH formation rely on stellar or IMBH mass seeds or direct collapse of gas clouds, and thus the detection or lack of IMBHs can help us understand SMBH formation (e.g., Greene et al. 2020).

Studying IMBHs and the lowest-mass SMBHs in galaxy centers can also help with extending and understanding the correlations that exist between galaxies and their black holes (e.g., Gebhardt et al. 2000; McConnell & Ma 2013; Saglia et al. 2016) to lower masses.

Recently, BHs with masses 10⁵–10⁷ $M_{\odot}$ have been detected in lower-mass galaxies with masses 10⁹–10¹⁰ $M_{\odot}$ using both dynamical measurements (den Brok et al. 2015; Nguyen et al. 2018, 2019; Davis et al. 2020), and measurements of AGNs (e.g., Reines et al. 2013; Chilingarian et al. 2018; Mezcua et al. 2018). SMBHs with masses >10⁶ $M_{\odot}$ have also been found at the centers of ultracompact dwarfs (UCDs; e.g., Seth et al. 2014; Ahn et al. 2017), massive star clusters that appear to be the tidally stripped nuclear star clusters of galaxies (e.g., Mieske et al. 2013; Pfeffer & Baumgardt 2013; Neumayer et al. 2020). While so far these BHs have only been found in the highest mass UCDs (Voggel et al. 2018), there are likely lower-mass stripped nuclei and BHs hiding among galaxies' globular cluster (GC) systems (Voggel et al. 2019). These objects are among the most likely targets for detecting IMBHs.

Although GCs are potential reservoirs for IMBHs, detecting these IMBHs remains challenging for several reasons. First, the gravitational sphere of influence of the IMBHs is small, which limits dynamical IMBH searches (that must resolve this radius) to within the Local Group. Second, dynamical evolution in GCs causes stellar-mass black holes (and more slowly, neutron stars) to mass segregate at the center of a cluster. Collections of these stellar remnants can create a rise in the central velocity dispersion, mimicking an IMBH (e.g., Baumgardt et al. 2019; Zocchi et al. 2019). While many stellar-mass BHs will be lost due to interactions or natal kicks, a significant fraction of BHs can be retained at the center in some clusters (up to ~2% of the cluster's mass; Weatherford et al. 2020). Lastly, radial anisotropy can also contribute to creating an observed rise in the central velocity dispersion without the presence of an IMBH (Zocchi et al. 2017).

Dynamical detections of IMBHs have been reported in the Milky Way, e.g., in ω Cen (Noyola et al. 2010; Baumgardt 2017).

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Liverpool John Moores University, UK</footnote>

<footnote n="2">²Department of Physics and Astronomy, University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA</footnote>

<footnote n="3">³Harvard-Smithsonian Center for Astrophysics, Observatory Building E, 60 Garden St., Cambridge, MA 02138, USA</footnote>

<footnote n="5">⁵Michigan State University, USA</footnote>

<footnote n="5">⁵Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany</footnote>

<footnote n="6">⁶European Space Agency, STScI, USA</footnote>

<footnote n="7">⁷Max Planck Institut für Astronomie, Heidelberg, Germany</footnote>

<footnote n="8">⁸Observatoire Astronomique de Strasbourg, USA</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Kinematics of B023-G78. The two panels are the stellar kinematic maps (velocity and velocity dispersion, respectively) of the cluster derived from adaptive-optics assisted Gemini/NIFS data. The systemic velocity ($V_{\text{sys}}$) was estimated to be -435 km s⁻¹.</figcaption></figure>

Due to the data's high spatial resolution, we were able to resolve individual stars in the GC. To mitigate shot-noise effects due to the brightest cluster stars, we used the *PampelMuse* software (Kamann et al. 2013) to generate a star-subtracted cube of the GC. To describe the method in short, we fitted a single Sérsic image to the continuum image of the cluster and inspected the fit residuals. We then manually identified the locations where the residuals suggested the presence of resolved stars and the resulting list was used to recover their PSFs as a function of wavelength using *PampelMuse* (Kamann et al. 2013). These PSF and the positions were used to extract the spectra of the input stellar sources. Finally, we combined the wavelength-dependent PSF model with the positions and the extracted spectra of the resolved stars to subtract their contributions from the NIFS data.

To derive the stellar kinematic maps, we first performed Voronoi binning using the code from Cappellari & Copin (2003). After resampling the integral-field spectroscopic data into bins of $S/N = 50$, we estimated the kinematics in each bin using the penalized pixel-fitting (pPXF) algorithm and code as described in Cappellari (2017). This code uses the full spectrum (1.5–1.8µ) to fit the radial velocity (V), and velocity dispersion ($\sigma_e$). We used 65 Phoenix stellar templates from Husser et al. (2013)<fnref n="10" /> with metallicities ranging from −1–0, log(g) from 1 to 5.5, temperatures from 3600 to 5500 K, and [α/Fe] from 0 to 0.4, covering the range of parameters expected to dominate the light in B023-G78. To estimate the kinematic errors, we performed Monte Carlo simulations by adding a random Gaussian error to the spectrum in each bin and re-fitting the kinematics. The standard deviation of those fits was taken as $1\sigma$ uncertainties. The final derived kinematics are shown in Figure 2. The central velocity dispersion is ~37 km s⁻¹, while clear rotation is seen around the minor axis with an amplitude of ~20 km s⁻¹ and a systemic velocity of ~435 km s⁻¹. The integrated dispersion out to 1" is 34.2 km s⁻¹; this is in reasonable agreement with the observed value of 31.7 ± 1.7 km s⁻¹ by Strader et al. (2011) using higher spectral resolution optical spectroscopy at seeing limited spatial resolution.

## 2.3. Deriving the Kinematic PSF

To perform dynamical modeling with precision, understanding the PSF of the Gemini/NIFS kinematic data is critical.

To determine the PSF, first we astrometrically aligned a continuum image created from the NIFS data cube (created without the additional on-chip sky subtraction) to the HST F814W image. Then, we used the HST F814W image within a 1" radius and convolved it with a double Gaussian model of the PSF, varying the parameters of the Gaussians until the convolved HST image best matched the continuum image created from the data cube. We then convolved the resulting PSF model with the HST PSF to obtain the widths and relative strengths of the two Gaussian components. The best-fit FWHM of the inner and outer Gaussian component was found to be 0.127 containing 31.4% of the total light, and 0.58 containing 68.6% of the light, respectively.

To account for systematic errors arising from the PSF model (described in Section 4.3), we also created PSFs with different sets of inputs. We estimated a PSF as described above, but fitted the F814W image out to a larger radius (1.35). Another PSF was generated as above but using the F606W HST image. The parameters of the PSFs in both the cases agreed within ~10% in the FWHMs and ratios of the components. This consistency is likely due to the lack of a strong color gradient seen in this cluster (see next section).

## 3. Creating Luminosity Models

We fit the HST-based surface brightness (SB) profiles of the cluster using the *IMFIT* code (Erwin 2015). The code builds a 2D model image using the input parameters and then convolves it with the given PSF. The best-fit $\chi^2$ is then estimated to find the closest model of the galaxy. Our best-fit model has two components, an inner King profile, and an outer Sérsic profile. The King profile (King 1962) is described by the central intensity ($I_0$), tidal radius ($r_t$), and core radius ($r_c$); although the *IMFIT* code uses a generalized King profile (Peng et al. 2010), we fix the power-law index $\alpha$ set to 2 to obtain the standard empirical King profile. The Sérsic profile (Sersic 1968; Graham & Driver 2005) is described by the Sérsic index (n), effective radius ($r_e$), and the intensity at the effective radius ($I_e$). Apart from these input parameters to the models, we also provide the position angle (PA), which is defined as the angle of the semimajor axis measured north through east counter-clockwise, and the ellipticity ($\epsilon$) as free parameters for each component to fit the 2D image of the cluster.

<footnote n="10">¹⁰ http://phoenix.astro.physik.uni-goettingen.de/</footnote>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Surface brightness profile and color of B023-G78. Left: best-fit model (King + Sérsic) decomposition of the surface brightness of B023-G78 in F814W and F606W from the ACS/HRC data; note that surface brightnesses are in Vega magnitudes. The blue (solid and dashed) is the M31 background estimated using SDSS images (in F814W and F606W, respectively; Section 3). The residuals (data—model) are shown in the bottom of the panel. Right: color map derived by convolving the PSF of F814W to the F606W image and vice versa (in black). The orange line is the color from the unconvolved model images, where the parameters of F814W were fixed to the best-fit model parameters of F606W, therefore requiring that the two components each have a unique color. The red line is the color when the fit parameters in both the F814W and F606W were allowed to be free. The magnitudes are extinction-corrected with an E(B-V) of 0.23.</figcaption></figure>

<table><caption>**Table 1**

Best-fit Parameters in F814W and F606W for B023-G78</caption><thead><tr><th>Function</th><th>Parameter</th><th>Best-fit Value F814W</th><th>Best-fit Values F606W</th></tr></thead><tbody><tr><td rowspan="5">King</td><td>log I<sub>0</sub></td><td>4.91 L<sub>☉</sub> pc<sup>-2</sup></td><td>4.65 L<sub>☉</sub> pc<sup>-2</sup></td></tr><tr><td>r<sub>c</sub></td><td>2.69 pc</td><td>2.68 pc</td></tr><tr><td>c = log(r<sub>f</sub>/r<sub>c</sub>)</td><td>1.11</td><td>1.12</td></tr><tr><td>ε</td><td>0.10</td><td>0.11</td></tr><tr><td>PA</td><td>80.0</td><td>76.4</td></tr><tr><td></td><td>mag<sub>tot</sub></td><td>13.02</td><td>14.22</td></tr><tr><td rowspan="6">Sérsic</td><td>log I<sub>e</sub></td><td>1.79 L<sub>☉</sub> pc<sup>-2</sup></td><td>1.24 L<sub>☉</sub> pc<sup>-2</sup></td></tr><tr><td>r<sub>e</sub></td><td>18.74 pc</td><td>15.06 pc</td></tr><tr><td>n</td><td>2.56</td><td>2.52</td></tr><tr><td>ε</td><td>0.24</td><td>0.26</td></tr><tr><td>PA</td><td>77.0</td><td>85.6</td></tr><tr><td>mag<sub>tot</sub></td><td>14.08</td><td>15.24</td></tr><tr><td>M31 background</td><td></td><td>21.57 mag arcsec<sup>-2</sup></td><td>22.59 mag arcsec<sup>-2</sup></td></tr><tr><td>Half-light radius</td><td>r<sub>hl</sub></td><td>4.23 pc</td><td></td></tr></tbody></table>

**Note.** The cluster is fitted by a King + Sérsic model. The parameters and their corresponding best-fit values are shown here. These parameters are from the free models. The magnitudes and luminosities are not extinction-corrected.

The image also includes background light from M31, which we assume is locally flat. We estimated the background from M31 and the sky levels of the HST images (which are not sky subtracted) by matching the SBs of the large-scale SDSS images in *r* and *i* bands transformed to HST F606W and F814W Vega magnitudes. After fitting for the HST image background levels, the transformed SDSS and HST SB profiles matched well beyond a radius of 5". To obtain the M31 background, we used the mean surface brightness value in the

region of 12"–18" from the SDSS image and incorporated this as a flat background component in our models. We estimated this background level to be at 21.57 and 22.59 mag arcsec<sup>-2</sup> in F606W and F814W, respectively. The standard deviation in the SDSS surface brightness at large radii was taken as the 1σ error on this estimate; ~0.05 mags in both bands.

The SB profile fits were performed separately in both the F814W and F606W filters. Here, the fitting parameters were allowed to vary. To estimate the change in the color of the cluster in the model images, we fixed the best-fit model parameters of the F814W image to that of the F606W image; we found an (F606W-F814W)<sub>0</sub> of 0.94 mags for the inner component and 0.89 mags for the outer component. A 1D radial profile of the cluster's surface brightness and the model fit is shown in the left panel of Figure 3, which was derived from summing up the 2D image in annuli of increasing radii. The best-fit parameters are given in Table 1. We note that Barmby et al. (2007) show that a single Wilson profile provides a good fit to the 1D profile of B023-G78. However, the combination of the outer component's bluer color and its significantly higher ellipticity indicates a real physical difference in the two components.

As noted earlier, the F814W PSF has a red halo, and thus to examine the color profile in more detail, we created a cross convolved color map, i.e., we convolved the F814W image with the PSF of F606W and vice versa. The resulting color profile is shown in the right panel of Figure 3. We show the central 5", out to the radius of our PSF. The observed color gradient roughly matches the expectations from our fixed parameter fits, with a ~0.05 mag decline between the central arcsecond and 5". Because this color gradient is so small, especially over the area we are fitting, we assume a constant *M/L* in our dynamical models, but we also explore mass models with a varying *M/L* in Section 4.3. The theoretical color for a population of 10 Gyr and [Fe/H] of -0.7 is ~0.82 mags using the PARSEC models (Bressan et al. 2012). This is considerably bluer than the observed and model colors,
<table><caption>**Table 2**
Best-fit MGE Parameters of F814W Fits for B023-G78</caption><thead><tr><th>Intensity<br>(L<sub>☉</sub> pc<sup>-2</sup>)</th><th>Gaussian Width<br>(arcsec)</th><th>Axial Ratio</th></tr></thead><tbody><tr><td>24526</td><td>0.19</td><td>0.90</td></tr><tr><td>2307</td><td>0.24</td><td>0.90</td></tr><tr><td>43610</td><td>0.40</td><td>0.90</td></tr><tr><td>11092</td><td>1.06</td><td>0.90</td></tr><tr><td>4396</td><td>0.07</td><td>0.76</td></tr><tr><td>3435</td><td>0.13</td><td>0.76</td></tr><tr><td>2278</td><td>0.23</td><td>0.76</td></tr><tr><td>694</td><td>0.29</td><td>0.76</td></tr><tr><td>1478</td><td>0.39</td><td>0.76</td></tr><tr><td>392</td><td>0.46</td><td>0.76</td></tr><tr><td>1169</td><td>0.62</td><td>0.76</td></tr><tr><td>805</td><td>0.90</td><td>0.76</td></tr><tr><td>252</td><td>1.14</td><td>0.76</td></tr><tr><td>190</td><td>1.34</td><td>0.76</td></tr><tr><td>454</td><td>1.67</td><td>0.76</td></tr><tr><td>288</td><td>2.80</td><td>0.76</td></tr><tr><td>1.40</td><td>3.13</td><td>0.76</td></tr><tr><td>84.4</td><td>5.77</td><td>0.76</td></tr></tbody></table>

Note. The MGE parameters for the fits to the F814W data in Table 1.

suggesting that our assumed reddening, $E(\text{B-V}) = 0.23$ (Jablonka et al. 1992), may be underestimated. We discuss this further in Section 4.2.

We use multi-Gaussian expansion (MGE) models to deproject the SB profiles for use in dynamical models. This method is described in detail in Pechetti et al. (2020). In short, we used the best-fit parameters from Table 2 and converted them to MGE models using the `mge_fit_1d` code (Cappellari 2002), sampling the SB profile logarithmically. The final model contains 18 Gaussian components as described in Table 2. The ellipticities from the `IMFIT` models were converted to axial ratios to deproject these MGEs.

## 4. Dynamical Modeling and BH Mass Estimates

In this section, we present dynamical models of B023-G078 that focus on constraining the mass of a possible central BH mass using Jeans' anisotropic modeling (JAM; Cappellari 2008). We first present results for our default model, then explore the impacts of the uncertain extinction correction and possible systematic errors on our best-fit models. We present additional dynamical models exploring the possibility of a cluster of stellar-mass black holes in Section 5.2.

### 4.1. Results from Jeans Anisotropic Modeling

For estimating the BH mass, we used the JAM method for our dynamical models. These models use the 3D deprojected MGE densities that were derived from the HST data in the previous section to create a gravitational potential. To this potential, a BH assuming a Gaussian potential with a very small scale (~0.001) is added. Using the potential and MGEs, the Jeans' equations are solved to estimate an intrinsic value of the root mean square (rms) velocity ($V_{\text{RMS}} = \sqrt{(V - V_{\text{sys}})^2 + \sigma_0^2}$), where V is the rotation velocity, $V_{\text{sys}}$ is the systemic velocity, and $\sigma_0$ is the velocity dispersion. The estimated $V_{\text{rms}}$ is then integrated along the line of sight to compare with the observed rms velocities derived from the Gemini/NIFS data out to a radius of 1". Our default model uses

the kinematic PSF derived from fitting the Gemini/NIFS data to the F814W image, the best-fit two-component King+Sérsic model derived from the F814W image (Table 1), and the kinematic data cube after star subtraction. We discuss additional models used to assess our systematic errors in Section 4.3.

We explore our JAM model fits by varying the following four free parameters: mass-to-light ratio ($M/L$), inclination angle (i), anisotropy parameter $\beta$, and BH mass $M_{\text{BH}}$, since they are degenerate. We estimate the best-fit values by sampling the parameter space using Markov Chain Monte Carlo (MCMC) simulations with the `emcee` package (Foreman-Mackey et al. 2013). We ran our models for 10,000 iterations. The resulting posterior probability distribution functions of our model parameters are shown in Figure 4.

We obtain a best-fit BH mass of $9.1_{-2.8}^{+2.6} \times 10^4 M_{\odot}$. The $\chi^2$ of the best-fit model is 404. The best-fit no-BH model has a $\Delta\chi^2$ of 30, excluding this model at $>3\sigma$ significance relative to the model with a BH. We estimated the Bayesian information criterion (BIC) for the best-fit IMBH model and the no-BH model. The $\Delta$BIC was 24, which provides strong evidence against the no-BH model. A $\Delta$BIC > 10 supports strong evidence for one model over another (Kass & Raftery 1995). For the best-fit BH mass and $\sigma_e$ as the integrated velocity dispersion at ~0.5, the sphere of influence radius (SOI = $GM_{\text{BH}}/\sigma_e^2$) is ~0.33 pc or ~0.09; for comparison, the PSF core sigma (FWHM/2.35) is 0.055, thus the SOI is resolved by our kinematic data as expected given the $>3\sigma$ significance of the BH mass detection. The best-fit $M/L_{\text{F814W}}$ is $1.87_{-0.04}^{+0.04}$, giving a total dynamical mass of $6.22 \times 10^6 M_{\odot}$ for this cluster. This total dynamical mass is similar to that found in Strader et al. (2011) ($6.8_{-0.6}^{+0.7} \times 10^6 M_{\odot}$). We discuss the uncertainties in the $M/L$ due to extinction in the next subsection but note that the dynamical mass is robust to changes in extinction. The models also suggest moderate radial anisotropy with a best-fit $\beta$ of $0.15_{-0.03}^{+0.05}$. The inclination is not well constrained, but does not affect the estimates of other parameters. To visualize the models and data better, Figure 5 shows a 1D radial profile of the measured annular $V_{\text{rms}}$ and the $V_{\text{rms}}$ model prediction, as well as showing the best-fit model without a BH. The model 1D profiles are estimated by creating radial bins from the 2D model and taking the median for each bin along the major axis of the cluster. Every iteration of the MCMC simulation within $1\sigma$ is also plotted, which is the shaded region.

### 4.2. Effects of Extinction

As noted in the 1, the dust lane passing in front of this cluster makes the extinction of this cluster poorly known. This uncertainty translates directly into an uncertainty in the $M/L$ of the cluster; however, the BH mass and total dynamical mass of the cluster are not sensitive to changes in the extinction because these are constrained by the combination of the kinematics and the shape of the mass model. To check for any extinction variations within the cluster, we averaged the color of the cluster azimuthally, but there were minimal variations (<1%). This suggests that the extinction is nearly constant in the region we are modeling. Our default reddening of $E(\text{B-V}) = 0.23$ (corresponding to $A_{\text{F814W}} = 0.427$) from Jablonka et al. (1992) provides an $M/L$ of ~1.9; the range of I band mass-to-light ratios observed in other M31 clusters is ~1-2 for high-mass clusters with similar metallicity (calculated from Peacock et al.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The output of JAM models from MCMC simulations showing the best-fit BH mass. $M_{BH}$ gives the black hole mass, $M/L$ indicates the mass-to-light ratio in the F814W band, $\beta$ shows the anisotropy, and $i$ gives the inclination. The top panel shows the probability distribution function of the black hole mass marginalized over all other parameters.</figcaption></figure>

2010; Strader et al. 2011), thus our derived value is reasonable, although a bit on the higher side. Using a higher reddening value in our dynamical modeling, like the $E(B-V) = 0.43$ from Caldwell et al. (2011), gives an $M/L_{F814W} = 1.35$, also still within the range of observed $M/Ls$. We also analyzed the resolved photometry in our HST data. Comparison of this data to Parsec isochrones (Bressan et al. 2012) suggests the CMD position of the RGB and red clump stars are consistent with $E(B-V) = 0.23$, and rules out significantly higher reddenings. We therefore use the Jablonka et al. (1992) value as our default value. We note that our choice of reddening/extinction does not impact the dynamical estimates of the best-fit BH mass or the total stellar mass of the cluster.

### 4.3. Sources of Systematic Error

Several systematic errors can affect our dynamical models. We discuss each of these and summarize their effect on our

estimated BH mass in the right panel of Figure 5. The default model as mentioned in the previous section is depicted in the black line.

One source of uncertainty in dynamically modeling GCs is defining the center (e.g., Noyola et al. 2010; Anderson & van der Marel 2010). When determining the surface photometry, *IMFIT* fits the center along with the cluster surface brightness. The formal error on the center was 0.156 mas, which underestimates the true uncertainty. We also estimated the center by running the ELLIPSE task using IRAF. We did not fix the center and estimated the center using ellipses with semimajor axes of $0.2''-3''$. We then determined the standard deviation of the measurements, which was ~12.5 mas. Given that this is ~1/4 the size of the kinematic pixels, this uncertainty has minimal impact on our dynamical models. Despite the small apparent uncertainty in our center, we tested the impact on the BH mass by shifting the central position of the cluster by 0.05
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: a 1D profile of the observed kinematics compared to model fits. Red points show annular $V_{rms}$ values from Gemini/NIFS. The black line shows the best-fit BH mass model with gray lines showing other models from the MCMC model fits. The blue line is the best-fit no-BH model derived from fitting only three parameters ($M/L$, *i*, and *β*). Right: cumulative likelihood of the BH mass estimate. The default model is the black solid line. The dashed lines are the $1\sigma$ uncertainty levels. We also show other models to highlight the level of systematic error in our default model. This includes mass-model variations (i.e., using the F606W image fits and varying the $M/L$ based on the color of the components), spatial offsets of one NIFS pixel, fits to kinematics derived from data cubes without star subtraction, and fits where we vary the kinematic PSF. All but the spatial offsets result in changes at the $<1\sigma$ level of the black hole mass, while the spatial offsets (which are significantly larger than our uncertainties in the center) are consistent within $2\sigma$ with the black hole mass in the default model.</figcaption></figure>

(1 NIFS pixel), in the *x* and *y* direction. The resulting variations in the cumulative distribution function of the inferred BH mass, shown as the red lines in Figure 5, were fairly large but still within the $2\sigma$ uncertainty of our default model (black line). Note that to get the center of kinematics to match that of the HST data, during our PSF analysis, we obtain the best-fit astrometry matching our NIFS data to the HST images.

Another major source of potential systematic error is the kinematic PSF that we derive using a double Gaussian profile. As described in Section 2.3, we estimated the PSFs using different fitting radii on the F814W image and using the F606W image. The impact of the PSF on our results is shown with green lines in Figure 5.

The luminosity/mass models we use also have two types of uncertainties: (1) uncertainties in the parameterization of the SB profiles, and (2) the possibility that the $M/L$ varies with radius, invalidating the mass-traces-light assumption in our first model. This could be due to varying stellar populations in the cluster (which we explore here) or due to mass segregation (discussed in the next subsection). To explore the size of uncertainties in (1) we use the best-fit F606W King+Sérsic model; this is shown as the blue line in Figure 5, which again did not create much variation in the BH mass. To explore (2), we performed tests by varying the $M/L$ in our mass model instead of assuming a constant $M/L$. We assigned a $M/L$ for the King and the Sérsic components based on their integrated colors using the theoretical color-$M/L$ relations from Roediger & Courteau (2015). These were then used as an input in our JAM models, and a single mass-scaling factor was used in place of the $M/L$ (as in Nguyen et al. 2018, 2019). This did not create much variation in the BH mass (shown as the pink line in Figure 5).

Finally, we used the original kinematics data cube rather than the one after star subtraction to estimate the BH mass but there was not much variation observed (shown as a cyan line in

Figure 5). Overall, these tests suggest that the dynamical signature of the IMBH in B023-G78 is robust.

Based on all the systematic errors that we explore, we find that none of them substantially change the estimated BH mass.

## 5. Discussion

We first discuss the evidence that B023-G078 is a stripped nucleus, and the interpretation of our BH results in that context. We then consider a collection of stellar-mass BHs as an alternative to the IMBH interpretation, and finish by examining B023-G078 in a broader context.

### 5.1. Additional Evidence That B023-G78 is a Stripped Nucleus

The presence of an IMBH might be expected in B023-G78 if it is a stripped nuclear star cluster (NSC) of a once more massive galaxy (e.g., Pfeffer & Baumgardt 2013). B023-G78 is the most massive cluster in M31 and an outlier in the M31 globular cluster luminosity function (Barmby et al. 2001; Strader et al. 2011). In the Milky Way, there is strong evidence that some of the most massive clusters are stripped NSCs. $\omega$Cen consists of complex stellar populations that cover a broad metallicity distribution. In addition, it has recently been suggested as the former core of Sequoia or the Gaia-Enceladus galaxy (e.g., Majewski et al. 2012; Myeong et al. 2019; Simpson et al. 2020; Pfeffer et al. 2021). The NSC of the tidally disrupting Sagittarius dwarf galaxy, M54, provides evidence that this stripping process is ongoing. It also shows complicated star formation history and a spread in metallicity and ages of the stars (e.g., Sarajedini & Layden 1995; Siegel et al. 2007; Alfaro-Cuello et al. 2019; Pfeffer et al. 2021). All of the globular cluster formation mechanisms are unable to explain these observations (e.g., Bastian & Lardo 2018). Assuming B023-G078 is a stripped NSC, we estimate a galaxy progenitor stellar mass of $5.3 \times 10^9 M_{\odot}$ using the galaxy-NSC
mass relation from Neumayer et al. (2020); this relation has
significant scatter, but most known NSCs of B023’s mass are
hosted in galaxies above 10⁹ M⊙. In this range of galaxy stellar
masses, the black hole occupation fraction is high (Nguyen
et al. 2019; Greene et al. 2020).

The metallicity (e.g., Janz et al. 2016) and metallicity spread (e.g., Pfeffer et al. 2021) of a globular cluster can provide additional evidence for a stripped NSC. The metallicity of B023-G078 has been estimated to be roughly -0.7 from several studies of both spectra and CMDs (Fuentes-Carrera et al. 2008; Perina et al. 2009; Caldwell et al. 2011). The observed metallicity for B023-G078 is within the observed range of NSC metallicities in its inferred host-galaxy mass range (Neumayer et al. 2020). Furthermore, the large spread in metallicities inferred from color–magnitude diagram modeling by Fuentes-Carrera et al. (2008) provides strong evidence for B023-G078 being a stripped nucleus with a range of metallicities similar to M54 (Alfaro-Cuello et al. 2019) or ω Cen (Johnson & Pilachowski 2010). The fits to APOGEE near-infrared spectra presented in Ashok et al. (2021) find a best-fit metallicity of -0.5±0.3, a best-fit [α/M] of +0.1, and a considerably younger age (~6 Gyr) than any of the other 32 M31 GCs analyzed; this younger population also is suggestive of a stripped NSC where young populations are expected to form until the epoch of stripping (Neumayer et al. 2020).

Given our observed color gradient and the previously
observed metallicity spread, we analyzed our NIFS spectra in
radial annuli to detect any significant age or metallicity gradient
in B023-G078. The spectra were binned into 12 annuli with a
maximum radius of 1"/25 with S/N ranging from >200 near
the center to ~80 at the largest radii. We then fit the spectra
using pPXF with the A-LIST spectral models, a set of simple
stellar population templates created using APOGEE spectra
(Ashok et al. 2021). We selected Padova-based templates with
[α/M] = 0.1, ages ranging from 2 to 12 Gyr, and [M/H] from
-2 to +0.4. The fits were very good, although due to the high
S/N, the reduced χ² was as high as 2.5 in the inner part of the
cluster. A light-weighted mean metallicity and age were
calculated at each radius and then a Monte Carlo analysis
was run to determine the errors on these quantities (note that
the error spectra were scaled by √χ² of the best-fit at each
radius during this analysis). The light-weighted metallicity is
consistent with previous metallicity determinations and shows
a clear negative gradient of ~0.15 dex between the center and
1" as shown in Figure 6. The light-weighted age is found to be
10.5 ± 0.5 Gyr with no significant gradient.<fnref n="11" /> The lower
metallicities at larger radii are also consistent with the bluer
colors of our outer component inferred in our model fits to the
B023-G078. The observed metallicity gradient is similar to
those seen in ω Cen and M54 (e.g., Suntzeff & Kraft 1996;
Monaco et al. 2005) with the metal-rich populations being
more concentrated than the metal-poor populations. Overall,
we interpret the metallicity spread and gradient as evidence of
the multiple generations of stars we expect to see in NSCs.

We note two additional pieces of evidence that favor B023-
G078 being a stripped NSC. First, the strong rotation
(V/σ≈0.8) seen is typical of NSCs (Neumayer et al. 2020),
but is higher than those seen in Milky Way GCs (Kamann et al.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** A clear metallicity gradient is seen in our Gemini/NIFS data. We fit annular spectra using A-LIST models (Ashok et al. 2021) to determine the light-weighted mean metallicity as a function of radius. Error bars are based on Monte Carlo simulations.</figcaption></figure>

2018); note that this value is a lower limit due to the unknown
inclination of the system. Second, the two-component structure
of the cluster is as expected from a stripped NSC (Pfeffer &
Baumgardt 2013), and is similar to the more massive UCDs
with known BHs (Seth et al. 2014; Ahn et al. 2017, 2018). The
apparent (weak) color variation between the two components is
also consistent with NSCs, where stellar population variations
and gradients are expected (Neumayer et al. 2020). Overall,
there is strong evidence that B023-G078 is in fact a stripped
nucleus from a galaxy in a mass range where central BHs are
commonly found.

5.2. Possible Alternatives to a Central IMBH

Dynamical evolution is expected to increase the *M/L* of
clusters both at the center, due to the mass segregation of BHs
and neutron stars, and in the outer parts, due to kicks received by
low-mass dwarf stars (e.g., den Brok et al. 2014; Baum-
gardt 2017). The mass segregation of the remnants happens on a
timescale less than the half-mass relaxation time, which in B023-
G78 is ~14 Gyr, and thus it is expected that the BH subsystem
will be mass-segregated. The expected mass fraction of stellar-
mass BHs retained over time remains extremely uncertain due to
poorly understood BH natal kicks from supernovae. Observa-
tionally, constraints on the BH kicks derived from the 3D
velocities of X-ray binaries suggest typical kicks >100 km s⁻¹
(Atri et al. 2019), with a small fraction having much lower kick
velocities; these are perhaps BHs formed from direct collapse.
The observed kicks are higher than expected from theoretical
prescriptions that base the natal kicks on the better constrained
neutron-star kick distribution with a linear decrease in mass due
to mass fall back and momentum conservation (e.g., Belczynski
et al. 2002; Morscher et al. 2015; Banerjee et al. 2020;
Mapelli 2021). The observed kick velocities are also above the
escape velocities of even the most massive Milky Way clusters,
including ω Cen (Gnedin et al. 2002). In addition to uncertainties
due to kicks, additional uncertainty on the retention fraction of
stellar-mass BHs comes from the unknown initial conditions for
clusters including the high-mass stellar initial mass function
(e.g., Baumgardt & Hilker 2018), and the uncertainty in the

<footnote n="11">¹¹ This age is significantly older than the value measured by Ashok et al. (2021) using the same models but independent data. We note that if we force a younger age on our fits, we get slightly worse fits, and higher metallicities consistent with the -0.5 found by Ashok et al. (2021).</footnote>
initial–final mass relation of BHs (e.g., Spera et al. 2015; Mapelli 2021).

Models with ~5% of the cluster mass in segregated stellar mass BHs are able to explain the rise in the central dispersion in ω Cen (Zocchi et al. 2019) and may be preferred to an IMBH model due to the lack of a high-velocity tail in the individual stellar velocities near the center (Baumgardt et al. 2019). An alternative constraint on BH mass fractions in Milky Way clusters was made by Weatherford et al. (2020) through modeling the observed mass segregation of stars. They constrained the BH mass fraction in Milky Way GCs, and found them to be <1% in 48/50 clusters (including the massive clusters 47 Tuc and M54). They report a correlation between the BH mass fraction and the ratio of the core radius to the half-light radius, and find two clusters with $r_c/r_{hl} > 0.75$ to have BH mass fractions of up to 2%. B023-G078’s large $r_c/r_{hl}$ thus suggests a high-mass fraction of BHs may be present. We note that the tidal stripping of star clusters can also lead to very high BH mass fractions as stars are lost from the cluster faster than the mass-segregated BHs (Gieles et al. 2021).

Relative to ω Cen, the higher metallicity of B023-G078 should lead to higher BH natal kicks (potentially lowering retention fractions) and lower typical BH masses and total BH mass fractions. To get a sense of the potential maximum mass fraction in BHs, we assumed a Kroupa (2001) IMF, the stellar evolution codes SSE & BSE using an [Fe/H] = −0.7 (Hurley et al. 2000, 2002), and the initial–final mass relations and BH kick prescriptions from Banerjee et al. (2020). This combination yields a total initial mass in BHs of 4.3%, making up 7.8% of the final mass. Removing BHs that receive kicks (and keeping only those that directly collapse), the present-day total mass fraction in BHs is 5.5%. As noted above, the retention of BHs is highly uncertain, and the kick prescription used here does not match that of observed X-ray binaries (Atri et al. 2019). However, the high mass of B023-G78 makes it plausible that a significant fraction of stellar-mass BHs are retained (e.g., Kremer et al. 2020). Thus it appears possible that the inferred central IMBH in B023-G78 may instead be a collection of stellar-mass BHs. We examine this possibility further below.

### 5.2.1. Testing the Stellar Mass BH Scenario with JAM Models

A collection of stellar-mass BHs differs from an IMBH because its mass distribution is extended, and this extent may be resolvable by our observations. In ω Cen, the best-fit distribution of stellar-mass BHs from Zocchi et al. (2019) can be described as a Plummer density profile ($\rho = 3M(1 + r^2/r_0^2)^{-5/2}/4\pi r_0^3$) with the ratio of the BH subsystem Plummer radius ($r_0$) to the cluster half-light radius ($r_{hl}$), ($r_0/r_{hl}$) ~ 0.3. This ratio would correspond to a $r_0$ of ~1.3 pc (0.′′3) in B023-G078; this is significantly broader than the core of our PSF and thus may result in measurable changes in our data relative to the point mass assumed in our IMBH models. However, we note that in the distribution function-based models of Zocchi et al. (2019), the amount of mass segregation between BHs and stars is fixed by a single parameter that is not well constrained, thus the ratio of ($r_0/r_{hl}$) is uncertain. A previous paper by Breen & Heggie (2013) used theory, gas models, and N-body models on idealized clusters to understand the expected distribution of their BHs; they found that for the parameters of ω Cen, the ratio of half-mass–radius of the BH subsystem ($r_{h,BH}$) over $r_{hl}$ is $r_{h,BH}/r_{hl} \simeq 0.15$. This ratio depends on the BH mass fraction; for a mass ratio of ~1% as we find for

the IMBH in B023-G078, they find $r_{h,BH}/r_h \simeq 0.1$. For a Plummer profile, this translates to $r_0/r_{hl} \sim 0.08$, which in B023-G078 would give an $r_0$ of just 0.3 pc, or 0.′′09, only slightly larger than the PSF core Gaussian width of 0.′′055 making for a more challenging measurement. Thus, if a significant BH subsystem is present in B023-G078, it is unclear whether we expect it to be significantly resolved by our observations.

To test whether an extended distribution of stellar mass BHs fits our kinematic data, we ran a new set of JAM models, replacing the central BH with a “dark” Plummer density profile. To include this in our JAM models, we created an MGE for the Plummer profile and included the Plummer radius ($r_0$) and the total mass (M) as free parameters in our MCMC simulations along with $M/L$, inclination, and $\beta$. The results are shown in Figure 7. From our simulations, we find that the median total mass of the dark component ($\sim 1.3 \times 10^5$) is within ~1$\sigma$ of our estimate for the IMBH. The median of the posterior for the Plummer $r_0$ parameter was ~0.11 pc (0.′′03) making it unresolved at our resolution; the best-fit value of 0.09 pc is also consistent with this small and unresolved $r_0$. The 95% confidence upper limit on $r_0$ is 0.56 pc, thus the upper limit on $r_0/r_{hl}$ is 0.13. The mass of the dark system increases with increasing size, and for the $r_0$ upper limit, the corresponding upper limit on the total mass of the BH subsystem is $2 \times 10^5 M_{\odot}$, 3.2% of the total system mass.

We also estimated the BIC for the IMBH simulations and the models with the Plummer profile. We find a ΔBIC of 6.3, providing positive evidence in favor of the IMBH models. Combining this with the considerable evidence that B023-G078 is a stripped NSC, we therefore favor an IMBH interpretation for our observations. However, a compact system of stellar-mass BHs is also a possible explanation for the observed rise in the central dispersion, as long as the $r_0 < 0.56$ pc and the total mass in the BH subsystem is $\lesssim 3\%$.

We note that it is possible that both an IMBH and a significant population of mass-segregated stellar-mass BHs are present. A central BH significantly slows the process of mass segregation but does not completely halt it (Antonini 2014). While we do not model this hybrid case here, the constraints on the total mass of the dark Plummer model above likely give an upper limit on the mass of the stellar mass BH subsystem, even in the case of co-existence with an IMBH. We summarize the results and the properties of B023-G078 in Tables 3 and 4.

### 5.3. B023-G78 in Context

Assuming our observed dynamical signature is an IMBH, we consider how it compares to other IMBH candidates and UCD/BH systems in Figure 8. At the lower-mass end, a comparison sample of claimed dynamical detections of massive BHs in GCs, as well as published upper limits for the same clusters, are shown from the recent compilation of Greene et al. (2020). We note many of the dynamical detections plotted here are disputed and refer readers to Greene et al. (2020) for details. In addition, we add higher-mass UCDs from recent discoveries, as well as present-day galaxies with both NSCs and BHs to provide context.

Relative to any other Local Group star cluster, the ~9 × $10^4 M_{\odot}$ BH in B023-G78 is the highest mass detection claimed, double the suggested mass of the BH in ω Cen (e.g., Noyola et al. 2010); as noted previously, this IMBH detection has been contested (e.g., Zocchi et al. 2017, 2019; Baumgardt et al. 2019). It is also more significant than the <3$\sigma$ detection
<figure><img src="image_10.png" /><figcaption>**Figure 7.** MCMC simulations of B023-G78 using a dark Plummer profile to describe a system of stellar mass BHs instead of an IMBH. The logM gives the total mass of the stellar mass BH subsystem, while logr0 indicates the Plummer radius; other parameters are the same as in Figure 4. The right two histograms give the best-fit “dark” component’s mass and size marginalized over all other parameters. The best-fit values of the total mass lies within 1σ of the IMBH mass in Figure 4, as do the inclination and anisotropy.</figcaption></figure>

of a $2 \times 10^4 M_{\odot}$ BH in G1 (Gebhardt et al. 2005) derived from data with similar physical resolution.

In comparison with the BHs previously found in other higher-mass UCDs, B023-G78 represents the first case in the IMBH regime, with all other BHs having both higher masses and mass fractions. Relative to central BHs in present-day galaxies, the mass is the lowest dynamical estimate apart from the ~$10^4 M_{\odot}$ BH suggested in NGC 205 (Nguyen et al. 2019). The most comparable present-day NSC+BH system is

NGC 4395, which hosts a ~$4 \times 10^5 M_{\odot}$ BH, inferred both dynamically (den Brok et al. 2015) and from reverberation mapping (e.g., Peterson et al. 2005), that lies in a ~$2 \times 10^6 M_{\odot}$ NSC (den Brok et al. 2015). The inferred IMBH in B023-G78 is also comparable to the lowest-mass BHs inferred from accretion (e.g., Baldassare et al. 2015; Chilingarian et al. 2018).

We also checked for possible BH accretion signatures in B023-G078. There is no cataloged X-ray source matching the location of B023-G078 in the deep XMM mosaic of
<table><caption>**Table 3**
Summary of Results</caption><thead><tr><th></th><th>IMBH model</th><th>No-BH model</th><th>Plummer model</th></tr></thead><tbody><tr><td><i>M</i><sub>BH</sub> (<i>M</i><sub>☉</sub>)</td><td>9.1<sup>+2.6</sup><sub>-2.8</sub> × 10<sup>4</sup></td><td>—</td><td>1.3<sup>+1.1</sup><sub>-0.4</sub> × 10<sup>5</sup></td></tr><tr><td><i>M</i>/L (<i>M</i><sub>☉</sub>/L<sub>☉</sub>)</td><td>1.87<sup>+0.04</sup><sub>-0.04</sub></td><td>1.99<sup>+0.02</sup><sub>-0.02</sub></td><td>1.83<sup>+0.05</sup><sub>-0.09</sub></td></tr><tr><td>β</td><td>0.15<sup>+0.06</sup><sub>-0.04</sub></td><td>0.22<sup>+0.04</sup><sub>-0.04</sub></td><td>0.16<sup>+0.04</sup><sub>-0.03</sub></td></tr><tr><td><i>i</i></td><td>64<sup>+18</sup><sub>-15</sub></td><td>59<sup>+16</sup><sub>-7</sub></td><td>64<sup>+13</sup><sub>-9</sub></td></tr><tr><td><i>r</i>0 (pc)</td><td>—</td><td>—</td><td>0.09</td></tr><tr><td><i>r</i>0 upper limit (pc)</td><td>...</td><td>...</td><td>0.56</td></tr><tr><td>Best-fit χ<sup>2</sup></td><td>404</td><td>434</td><td>404</td></tr></tbody></table>

**Note.** The best-fit parameters from the three different models we fit to the cluster.

<table><caption>**Table 4**
B023-G078 Cluster Properties</caption><tr><td>Central <i>V</i><sub>rms</sub></td><td>37.2 ± 0.6 km s<sup>-1</sup></td></tr><tr><td><i>V</i>/<i>σ</i></td><td>0.8</td></tr><tr><td>Cluster mass</td><td>6.22<sup>+0.03</sup><sub>-0.02</sub> × 10<sup>6</sup> M<sub>☉</sub></td></tr><tr><td>BH Mass</td><td>9.1<sup>+2.6</sup><sub>-2.8</sub> × 10<sup>4</sup> M<sub>☉</sub></td></tr><tr><td>BH mass fraction</td><td>1.5%</td></tr><tr><td>Half-mass relaxation time</td><td>14 Gyr</td></tr><tr><td>[Fe/H]</td><td>-0.65 (center) to -0.80 (at 1")</td></tr><tr><td>Age</td><td>10.5 ± 0.5 Gyr</td></tr><tr><td>Assumed E(B-V)</td><td>0.23</td></tr></table>

**Note.** B023-G078 properties that we derived from our analyses. The E(B-V) value is from Jablonka et al. (1992) and used as a default value in this paper.

Stiele et al. (2011). The faintest cataloged sources close to the location of B023-G078 have 0.5–4.5 keV XMM/EPIC unabsorbed fluxes of $2.1 \times 10^{-14}$ erg s$^{-1}$ cm$^{-2}$, which corresponds to a 0.5–10 keV X-ray luminosity of $1.9 \times 10^{35}$ erg$^{-1}$ assuming a photon index of $\Gamma = 1.7$. Hence the nondetection of B023-G078 in these data suggests a 0.5–10 keV upper limit of $L_X \lesssim 2 \times 10^{35}$ erg$^{-1}$. Using this upper limit to the X-ray luminosity in B023-G078 combined with our derived dynamical BH mass, the predicted 5 GHz luminosity is $< 8.5 \mu Jy$ (Plotkin et al. 2012). B023-G078 is not detected in VLASS, and with an rms noise in the VLASS image of 127 $\mu Jy$/bm, we can estimate a $3\sigma$ upper limit of $< 381 \mu Jy$. Therefore, in this case, the X-ray limit (if accurate) is much more constraining than the radio data, although it would be possible to get significantly deeper radio data. We also note that the presence of stellar mass black holes could also lead to detectable X-ray binaries, as B023-G078 does have a very high collision rate. However, among the highest collision rate GCs in M31 only a fraction ($< \text{half}$) appear to have bright X-ray sources (e.g., Peacock et al. 2010). We note in this context that in $\omega$ Cen, which as discussed above, may host a large cluster of stellar mass BHs (Zocchi et al. 2019; Baumgardt et al. 2019). However, no bright X-ray binaries are found, with the brightest X-ray sources being $< 10^{33}$ erg s$^{-1}$ (Henleywillis et al. 2018).

One potentially comparable system detected via accretion is HLX-1, a bright off-nuclear X-ray source with an inferred BH mass of $10^{4-5} M_{\odot}$ (e.g., Webb et al. 2012). Due to the light from HLX-1 itself, constraining the age and mass of the surrounding stellar cluster is challenging (e.g., Soria et al. 2010; Farrell et al. 2014; Soria et al. 2017), but if it is old, its mass is estimated to be $\sim 3 \times 10^6 M_{\odot}$ (Soria et al. 2017).

<figure><img src="image_11.png" /><figcaption>**Figure 8.** The black hole—cluster mass diagram for systems comparable to B023-G078BH. The blue data points are BH mass estimates in GCs from the compilation of Greene et al. (2020). The green data points are the fundamental-plane upper limits of the same clusters from the same compilation. The orange points are stripped nuclei (UCDs) from Seth et al. (2014), Ahn et al. (2017, 2018), Voggel et al. (2018), and Afanasiev et al. (2018). The gray points are the objects that have an estimate of both the BH and NSC mass from the compilation in Neumayer et al. (2020). All the open circles are upper limits. The dashed lines show BH masses that are 1% and 10% of their cluster mass. B023-G78 is the highest mass BH detected in a cluster below $10^7 M_{\odot}$.</figcaption></figure>

## 6. Conclusions

We have presented adaptive-optics GEMINI/NIFS IFU kinematic data of M31’s most massive star cluster, B023-G78. We combined these data with mass models derived from HST ACS/HRC to constrain the mass content, including a possible central black hole in this massive star cluster. We find the following:

1. The kinematics of B023-G78 show a rise in the integrated velocity dispersion to $\sim 37$ km s$^{-1}$, and a peak rotation of $\sim 20$ km s$^{-1}$.

2. The surface brightness profile requires at least two components to fit, and shows a small color gradient, with the outer component being $\sim 0.05$ mag bluer than the inner component. A significant metallicity gradient of $\sim 0.15$ dex is also seen within the central arcsecond.

3. Our best-fit JAM dynamical models give a BH mass of $9.1_{-2.8}^{+2.6} \times 10^4 M_{\odot}$, $M/L_{F814W}$ of $1.87_{-0.04}^{+0.04}$ and anisotropy $0.15_{-0.04}^{+0.06}$. The BH detection is highly significant $> 3\sigma$, and systematic errors are $< 10\%$ on the best-fit BH mass.

We discuss the possibility that this BH can be explained due to a collection of dark stellar remnants and constrain the extent of these remnants and find the derived extent of the dark remnants are mostly unresolved by our observations, with an upper limit on the Plummer $r_0$ of 0.56 pc. We favor the presence of a single IMBH given the other indications that B023-G78 is a stripped nucleus, as well as the apparent compactness of the dark component. Higher spatial-resolution data would give improved constraints on the nature of the central dark mass and should be a high priority in the forthcoming era of extremely large telescopes.
The authors thank Mark Gieles and Alice Zocchi for useful conversations about this work. R.P. and A.C.S acknowledge support from grants NSF AST-1350389 and AST-1813708. N.N. gratefully acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 138713538—SFB 881 (“The Milky Way System,” subproject B8). R.P. and S.K. gratefully acknowledge funding from UKRI in the form of a Future Leaders Fellowship (grant no. MR/T022868/1. J.S. acknowledges support from NSF grant AST-1812856 and the Packard Foundation. This paper is based on observations obtained at the international Gemini Observatory, a program of NSFs NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).

## ORCID iDs

Renuka Pechetti https://orcid.org/0000-0002-1670-0808

Anil Seth https://orcid.org/0000-0003-0248-5470

Sebastian Kamann https://orcid.org/0000-0001-6604-0505

Nelson Caldwell https://orcid.org/0000-0003-2352-3202

Jay Strader https://orcid.org/0000-0002-1468-9668

Nora Luetzgendorf https://orcid.org/0000-0002-4034-0080

Nadine Neumayer https://orcid.org/0000-0002-6922-2598

Karina Voggel https://orcid.org/0000-0001-6215-0950

## References

* Afanasiev, A. V., Chilingarian, I. V., Mieske, S., et al. 2018, MNRAS, 477, 4856

* Ahn, C. P., Seth, A. C., Cappellari, M., et al. 2018, ApJ, 858, 102

* Ahn, C. P., Seth, A. C., den Brok, M., et al. 2017, ApJ, 839, 72

* Alfaro-Cuello, M., Kacharov, N., Neumayer, N., et al. 2019, ApJ, 886, 57

* Anderson, J., & van der Marel, R. P. 2010, ApJ, 710, 1032

* Antonini, F. 2014, ApJ, 794, 106

* Ashok, A., Zasowski, G., Seth, A., et al. 2021, AJ, 161, 167

* Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2019, MNRAS, 489, 3116

* Baldassare, V. F., Reines, A. E., Gallo, E., & Greene, J. E. 2015, ApJL, 809, L14

* Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41

* Barnby, P., Huchra, J. P., & Brodie, J. P. 2001, AJ, 121, 1482

* Barnby, P., McLaughlin, D. E., Harris, W. E., Harris, G. L. H., & Forbes, D. A. 2007, AJ, 133, 2764

* Bastian, N., & Lardo, C. 2018, ARA&A, 56, 83

* Baumgardt, H. 2017, MNRAS, 464, 2174

* Baumgardt, H., He, C., Sweet, S. M., et al. 2019, MNRAS, 488, 5340

* Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520

* Baumgardt, H., Makino, J., Hut, P., McMillan, S., & Portegies Zwart, S. 2003, ApJL, 589, L25

* Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407

* Breen, P. G., & Heggie, D. C. 2013, MNRAS, 432, 2779

* Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

* Caldwell, N., Schiavon, R., Morrison, H., Rose, J. A., & Harding, P. 2011, AJ, 141, 61

* Cappellari, M. 2002, MNRAS, 333, 400

* Cappellari, M. 2008, MNRAS, 390, 71

* Cappellari, M. 2017, MNRAS, 466, 798

* Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345

* Chilingarian, I. V., Katkov, I. Y., Zolotukhin, I. Y., et al. 2018, ApJ, 863, 1

* Davis, S. W., Narayan, R., Zhu, Y., et al. 2011, ApJ, 734, 111

* Davis, T. A., Nguyen, D. D., Seth, A. C., et al. 2020, MNRAS, 496, 4061

* den Brok, M., Seth, A. C., Barth, A. J., et al. 2015, ApJ, 809, 101

* den Brok, M., van de Ven, G., van den Bosch, R., & Watkins, L. 2014, MNRAS, 438, 487

* Erwin, P. 2015, ApJ, 799, 226

* Farrell, S. A., Servillat, M., Gladstone, J. C., et al. 2014, MNRAS, 437, 1208

* Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

* Fuentes-Carrera, I., Jablonka, P., Sarajedini, A., et al. 2008, A&A, 483, 769

* Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJL, 539, L13

* Gebhardt, K., Rich, R. M., & Ho, L. C. 2002, ApJL, 578, L41

* Gebhardt, K., Rich, R. M., & Ho, L. C. 2005, ApJ, 634, 1093

* Gieles, M., Erkal, D., Antonini, F., Balbinot, E., & Peñarrubia, J. 2021, NatAs, 5, 957

* Gnedin, O. Y., Zhao, H., Pringle, J. E., et al. 2002, ApJL, 568, L23

* Godet, O., Plazolles, B., Kawaguchi, T., et al. 2012, ApJ, 752, 34

* Gonzaga, S., et al. 2012, The DrizzlePac Handbook (Baltimore, MD: STScI)

* Graham, A. W., & Driver, S. P. 2005, PASA, 22, 118

* Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257

* Henleywillis, S., Cool, A. M., Haggard, D., et al. 2018, MNRAS, 479, 2834

* Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543

* Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897

* Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

* Ibata, R., Bellazzini, M., Chapman, S. C., et al. 2009, ApJL, 699, L169

* Jablonka, P., Alloin, D., & Bica, E. 1992, A&A, 260, 97

* Janz, J., Norris, M. A., Forbes, D. A., et al. 2016, MNRAS, 456, 617

* Johnson, C. I., & Pilachowski, C. A. 2010, ApJ, 722, 1373

* Kamann, S., Husser, T. O., & Dreizler, S. 2018, MNRAS, 473, 5591

* Kamann, S., Wisotzki, L., & Roth, M. M. 2013, A&A, 549, A71

* Karachentsev, I. D., Karachentseva, V. E., Huchtmeier, W. K., & Makarov, D. I. 2004, AJ, 127, 2031

* Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

* King, I. 1962, AJ, 67, 471

* Kremer, K., Ye, C. S., Rui, N. Z., et al. 2020, ApJS, 247, 48

* Kroupa, P. 2001, MNRAS, 322, 231

* Lützgendorf, N., Gebhardt, K., Baumgardt, H., et al. 2015, A&A, 581, A1

* Majewski, S. R., Nidever, D. L., Smith, V. V., et al. 2012, ApJL, 747, L37

* Mapelli, M. 2021, in Formation Channels of Single and Binary Stellar-Mass Black Holes, ed. C. Bambi, S. Katsanevas, & K. D. Kokkotas (Singapore: Springer)

* McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184

* Mezcua, M., Civano, F., Marchesi, S., et al. 2018, MNRAS, 478, 2576

* Mieske, S., Frank, M. J., Baumgardt, H., et al. 2013, A&A, 558, A14

* Miller-Jones, J. C. A., Wrobel, J. M., Sivakoff, G. R., et al. 2012, ApJL, 755, L1

* Monaco, L., Bellazzini, M., Ferraro, F. R., & Pancino, E. 2005, MNRAS, 356, 1396

* Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9

* Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235

* Neumayer, N., Seth, A., & Böker, T. 2020, A&Ar, 28, 4

* Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 118

* Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2019, ApJ, 872, 104

* Noyola, E., Gebhardt, K., Kissler-Patig, M., et al. 2010, ApJL, 719, L60

* Peacock, M. B., Maccarone, T. J., Knigge, C., et al. 2010, MNRAS, 402, 803

* Pechetti, R., Seth, A., Neumayer, N., et al. 2020, ApJ, 900, 32

* Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

* Perina, S., Federici, L., Bellazzini, M., et al. 2009, A&A, 507, 1375

* Peterson, B. M., Bentz, M. C., Desroches, L.-B., et al. 2005, ApJ, 632, 799

* Pfeffer, J., & Baumgardt, H. 2013, MNRAS, 433, 1997

* Pfeffer, J., Lardo, C., Bastian, N., Saracino, S., & Kamann, S. 2021, MNRAS, 500, 2514

* Plotkin, R. M., Markoff, S., Kelly, B. C., Körding, E., & Anderson, S. F. 2012, MNRAS, 419, 267

* Reines, A. E., Greene, J. E., & Geha, M. 2013, ApJ, 775, 116

* Roediger, J. C., & Courteau, S. 2015, MNRAS, 452, 3209

* Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47

* Sarajedini, A., & Layden, A. C. 1995, AJ, 109, 1086

* Sersic, J. L. 1968, Atlas de Galaxias Australes (Cordoba: Observatorio Astronomico)

* Seth, A. C., Cappellari, M., Neumayer, N., et al. 2010, ApJ, 714, 713

* Seth, A. C., van den Bosch, R., Mieske, S., et al. 2014, Natur, 513, 398

* Siegel, M. H., Dotter, A., Majewski, S. R., et al. 2007, ApJL, 667, L57

* Simpson, J. D., Martell, S. L., Da Costa, G., et al. 2020, MNRAS, 491, 3374

* Soria, R., Hau, G. K. T., Graham, A. W., et al. 2010, MNRAS, 405, 870
Soria, R., Musaeva, A., Wu, K., et al. 2017, MNRAS, 469, 886

Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086

Stiele, H., Pietsch, W., Haberl, F., et al. 2011, A&A, 534, A55

Strader, J., Caldwell, N., & Seth, A. C. 2011, AJ, 142, 8

Straub, O., Godet, O., Webb, N., Servillat, M., & Barret, D. 2014, A&A, 569, A116

Suntzeff, N. B., & Kraft, R. P. 1996, AJ, 111, 1913

Tremou, E., Strader, J., Chomiuk, L., et al. 2018, ApJ, 862, 16

Voggel, K. T., Seth, A. C., Baumgardt, H., et al. 2019, ApJ, 871, 159

Voggel, K. T., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 20

Weatherford, N. C., Chatterjee, S., Kremer, K., & Rasio, F. A. 2020, ApJ, 898, 162

Webb, N., Cseh, D., Lenc, E., et al. 2012, Sci, 337, 554

Zocchi, A., Gieles, M., & Hénault-Brunet, V. 2017, MNRAS, 468, 4429

Zocchi, A., Gieles, M., & Hénault-Brunet, V. 2019, MNRAS, 482, 4713