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

We used archival HST data for this cluster from the proposal ID:9719 (PI: T. Bridges). ⁹ The observations were performed with the Advanced Camera for Surveys/High-Resolution Camera (ACS/HRC) in the filters F814W and F606W. The exposure times were 2860 s and 2020 s, respectively.

The ACS/HRC has a pixel scale of 0025 pixel⁻¹ and a field of view of 29'' × 26''. We downloaded the individual.flt files from the Mikulski archive for space telescopes (MAST) and drizzled them using Astrodrizzle (Gonzaga et al. 2012) to create the final image. The MDRIZSKY keyword was set to zero to avoid oversubtraction of the sky in the final drizzled image. The final color image (F606W—F814W) is shown in the right panel of Figure 1.

Zoom In Zoom Out Reset image size

The drizzled PSF in each band was obtained by inserting Tiny Tim PSFs into mock.flt images and drizzling them the same way as the science image. This procedure is similar to the one described in Pechetti et al. (2020). We note that the F814W PSF has a significant amount of light in a large halo; a PSF of radius 5'' was used to account for this.

We obtained Gemini/NIFS laser-guide-star adaptive-optics observations of B023-G78 on 2014 October 7 and November 9 as part of program GN-2014B-DD-2 (PI: A.C. Seth). The data provide integral-field spectroscopy in the H band (1.48–1.79 μm) over a field of view of 3'' (11 pc at 0.77 Mpc). For our final data cube with a spaxel size of 005, we combined the 6/8 900 s dithered exposures using the Gemini IRAF packages, with modifications as described in Seth et al. (2010) and Ahn et al. (2018). Despite the use of offset sky exposures, additional on-chip sky subtraction was required before combination, using the corners of the chip; this makes our useful field of view ∼2''. The line spread function was measured from sky lines in each pixel with a median FWHM of 3.27 Å. The kinematic maps are shown in Figure 2.

Zoom In Zoom Out Reset image size

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 (σ_e ). We used 65 Phoenix stellar templates from Husser et al. (2013) ¹⁰ 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σ 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.

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 0127 containing 31.4% of the total light, and 058 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 (135). 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).

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 (∼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}_{{RMS}}=\sqrt{{(V-{V}_{\mathrm{sys}})}^{2}+{\sigma }_{0}^{2}}$), where V is the rotation velocity, V_sys is the systemic velocity, and σ₀ is the velocity dispersion. The estimated V_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 β, and BH mass M_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.

Zoom In Zoom Out Reset image size

We obtain a best-fit BH mass of ${9.1}_{-2.8}^{+2.6}\times {10}^{4}$ M_⊙. The χ² of the best-fit model is 404. The best-fit no-BH model has a Δχ² of 30, excluding this model at >3σ 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 ΔBIC was 24, which provides strong evidence against the no-BH model. A ΔBIC > 10 supports strong evidence for one model over another (Kass & Raftery 1995). For the best-fit BH mass and σ_e as the integrated velocity dispersion at ∼05, the sphere of influence radius ($\mathrm{SOI}={{GM}}_{{BH}}/{\sigma }_{e}^{2}$) is ∼0.33 pc or ∼009; 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σ significance of the BH mass detection. The best-fit M/L_F814W is ${1.87}_{-0.04}^{+0.04}$, giving a total dynamical mass of 6.22 × 10⁶ M_⊙ 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_⊙). 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 β 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_rms and the V_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σ is also plotted, which is the shaded region.

Zoom In Zoom Out Reset image size

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(B–V) = 0.23 (corresponding to A_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. 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.

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 02–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 005 (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σ 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.

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. ω 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 × 10⁹ M_⊙ 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.1}^{+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 125 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 $\sqrt{{\chi }^{2}}$ 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. ¹¹ 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.

Zoom In Zoom Out Reset image size

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. 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.

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; Baumgardt 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. Observationally, 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 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₀) to the cluster half-light radius (r_hl), (r₀/r_hl) ∼ 0.3. This ratio would correspond to a r₀ of ∼1.3 pc (03) 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₀/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 ≃ 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 ≃ 0.1. For a Plummer profile, this translates to r₀/r_hl ∼ 0.08, which in B023-G078 would give an r₀ of just 0.3 pc, or 009, only slightly larger than the PSF core Gaussian width of 0055 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₀) and the total mass (M) as free parameters in our MCMC simulations along with M/L, inclination, and β. The results are shown in Figure 7. From our simulations, we find that the median total mass of the dark component (∼1.3 × 10⁵) is within ∼1σ of our estimate for the IMBH. The median of the posterior for the Plummer r₀ parameter was ∼0.11 pc (003) making it unresolved at our resolution; the best-fit value of 0.09 pc is also consistent with this small and unresolved r₀. The 95% confidence upper limit on r₀ is 0.56 pc, thus the upper limit on r₀/r_hl is 0.13. The mass of the dark system increases with increasing size, and for the r₀ upper limit, the corresponding upper limit on the total mass of the BH subsystem is 2 × 10⁵ M_⊙, 3.2% of the total system mass.

Zoom In Zoom Out Reset image size

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.56 pc and the total mass in the BH subsystem is ≲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.

Table 3. Summary of Results

	IMBH	No-BH	Plummer
	model	model	model
MBH (M⊙)	${9.1}_{-2.8}^{+2.6}\times {10}^{4}$	−	${1.3}_{-0.4}^{+1.1}\times {10}^{5}$
M/L (M⊙/L⊙)	${1.87}_{-0.04}^{+0.04}$	${1.99}_{-0.02}^{+0.02}$	1.83${}_{-0.09}^{+0.05}$
β	${0.15}_{-0.04}^{+0.06}$	${0.22}_{-0.04}^{+0.04}$	${0.16}_{-0.03}^{+0.04}$
i	${64}_{-15}^{+18}$	${59}_{-7}^{+16}$	${64}_{-9}^{+13}$
r0 (pc)	−	−	0.09
r0 upper limit (pc)	⋯	⋯	0.56
Best-fit χ2	404	434	404

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

Download table as:  ASCIITypeset image

Table 4. B023-G078 Cluster Properties

Central Vrms	37.2 ± 0.6 km s−1
V/σ	0.8
Cluster mass	6.22${}_{-0.02}^{+0.03}\times {10}^{6}$ M⊙
BH Mass	${9.1}_{-2.8}^{+2.6}\times {10}^{4}$ M⊙
BH mass fraction	1.5%
Half-mass relaxation time	14 Gyr
[Fe/H]	−0.65 (center) to −0.80 (at 1'')
Age	10.5 ± 0.5 Gyr
Assumed E(B–V)	0.23

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.

Download table as:  ASCIITypeset image

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.

Zoom In Zoom Out Reset image size

Relative to any other Local Group star cluster, the ∼9 × 10⁴ M_⊙ 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σ detection of a 2 × 10⁴ M_⊙ 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⁴ M_⊙ BH suggested in NGC 205 (Nguyen et al. 2019). The most comparable present-day NSC+BH system is NGC 4395, which hosts a ∼4 × 10⁵ M_⊙ BH, inferred both dynamically (den Brok et al. 2015) and from reverberation mapping (e.g., Peterson et al. 2005), that lies in a ∼2 × 10⁶ M_⊙ 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 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 × 10⁻¹⁴ erg s⁻¹ cm⁻², which corresponds to a 0.5–10 keV X-ray luminosity of 1.9 × 10³⁵ erg ⁻¹ assuming a photon index of Γ = 1.7. Hence the nondetection of B023-G078 in these data suggests a 0.5–10 keV upper limit of L_X ≲ 2 × 10³⁵ erg ⁻¹. 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 μJy (Plotkin et al. 2012). B023-G078 is not detected in VLASS, and with an rms noise in the VLASS image of 127 μJy/bm, we can estimate a 3σ upper limit of < 381 μ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 (< half) appear to have bright X-ray sources (e.g., Peacock et al. 2010). We note in this context that in ω 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³³ erg s⁻¹ (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_⊙ (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 ∼3 × 10⁶ M_⊙ (Soria et al. 2017).