A 40-billion solar mass black hole in the extreme core of Holm 15A, the central galaxy of Abell 85

TL;DR

This study reveals Holm 15A hosts the most massive black hole ever directly measured in the local universe, using advanced dynamical modeling of stellar kinematics, and links its properties to galaxy merger history.

Contribution

It provides the first direct dynamical measurement of a supermassive black hole in Holm 15A, the largest in the local universe, and explores its formation through galaxy mergers.

Findings

Holm 15A's black hole mass is (4.0 +- 0.80) x 10^10 solar masses.

The galaxy shows a tangential orbital bias less than other cored ellipticals.

Black hole mass correlates inversely with central stellar surface brightness and density.

Abstract

Holm 15A, the brightest cluster galaxy (BCG) of the galaxy cluster Abell 85, has an ultra-diffuse central region, 2 mag fainter than the faintest depleted core of any early-type galaxy (ETG) that has been dynamically modelled in detail. We use orbit-based, axisymmetric Schwarzschild models to analyse the stellar kinematics of Holm 15A from new high-resolution, wide-field spectral observations obtained with MUSE at the VLT. We find a supermassive black hole (SMBH) with a mass of (4.0 +- 0.80) x 10^10 solar masses at the center of Holm 15A. This is the most massive black hole with a direct dynamical detection in the local universe. We find that the distribution of stellar orbits is increasingly biased towards tangential motions inside the core. However, the tangential bias is less than in other cored elliptical galaxies. We compare Holm 15A with N-body simulations of mergers between galaxies with black holes and find that the observed amount of tangential anisotropy and the shape of the light profile are consistent with a formation scenario where Holm 15A is the remnant of a merger between two ETGs with pre-existing depleted cores. We find that black hole masses in cored galaxies, including Holm 15A, scale inversely with the central stellar surface brightness and mass density, respectively. These correlation are independent of a specific parameterization of the light profile.

Figures

Tables

Table 1. Table 1: Results of Schwarzschild dynamical modeling of Holm 15A. Best-fit values were derived as the mean of the independent fits to the four quadrants. The quoted uncertainties are derived from the variation between quadrants.

Schwarzschild Model Parameter	Best-Fit Value	Units
\toprule$M_{BH}$	$\left(4.0\pm 0.80\right)$	${10}^{10}{M}_{\odot }$
${{\rm Y}}_{\star }$ (i-band)	$4.5\pm 0.19$
DM Halo:
${\rho }_{10}$	$\left(1.0\pm 0.10\right)$	${10}^7\frac{M_{\odot }}{kp{c}^3}$
$\log r_s$	$\left(2.4\pm 0.29\right)$	$\log \frac{r}{kpc}$
$\gamma$	$0.35\pm 0.26$

Table 2. Table 2: Parameters of our best-fit models to the 1D i 𝑖 i -band Wendelstein image of Holm 15A, separated into components: cSS: Core-Sérsic + (outer) Sérsic fit to the light rofile out to 200   ″ times 200 ″ 200\text{\,}\arcsec with all parameters fit simultaneously. cSS( r m ​ i ​ n = 4 ​ ″ subscript 𝑟 𝑚 𝑖 𝑛 4 ″ r_{min}=4\arcsec ): Core-Sérsic + (outer) Sérsic fit to the light profile out to 200   ″ times 200 ″ 200\text{\,}\arcsec but with the parameters fit in two steps as described in Section 2.2 with r m ​ i ​ n = 4 ​ ″ subscript 𝑟 𝑚 𝑖 𝑛 4 ″ r_{min}=4\arcsec . cSS( r m ​ i ​ n = 12 ​ ″ subscript 𝑟 𝑚 𝑖 𝑛 12 ″ r_{min}=12\arcsec ): same as the previous model, but with r m ​ i ​ n = 12 ​ ″ subscript 𝑟 𝑚 𝑖 𝑛 12 ″ r_{min}=12\arcsec . N( r m ​ a ​ x = 20 ​ ″ subscript 𝑟 𝑚 𝑎 𝑥 20 ″ r_{max}=20\arcsec ): Nuker profile fit to the the data within r m ​ a ​ x = 20 ​ ″ subscript 𝑟 𝑚 𝑎 𝑥 20 ″ r_{max}=20\arcsec with all 5 parameters simultaneously. N( r m ​ a ​ x = 70 ​ ″ subscript 𝑟 𝑚 𝑎 𝑥 70 ″ r_{max}=70\arcsec ): same as the previous model, but with r m ​ a ​ x = 70 ​ ″ subscript 𝑟 𝑚 𝑎 𝑥 70 ″ r_{max}=70\arcsec .

Model	Parameter	cSS	cSS($r_{min}=4″$)	cSS($r_{min}=12″$)	N($r_{max}=20″$)	N($r_{max}=70″$)	Units
Core-Sérsic	$r_b$	$8.96$	$3.72$	$5.52$	$\dots$	$\dots$	arcsec
	${\mu }_b$	$20.6$	$20.1$	$20.1$	$\dots$	$\dots$	mag arcsec-2
	$\alpha$	$1.71$	$7.96$	$2.65$	$\dots$	$\dots$
	$\gamma$	$0.09$	$0.00$	$0.13$	$\dots$	$\dots$
	$n_1$	$5.24$	$1.38$	$2.52$	$\dots$	$\dots$
	$r_{e,1}$	$5.20$	$14.4$	$14.9$
Outer Sérsic	$n_2$	$2.90$	$5.37$	$3.30$	$\dots$	$\dots$
	${\mu }_{e,2}$	$28.7$	$28.6$	$28.8$	$\dots$	$\dots$	mag arcsec-2
	$r_{e,2}$	$917.0$	$647.3$	$952.5$	$\dots$	$\dots$	arcsec
Nuker	$r_b$	$\dots$	$\dots$	$\dots$	$11.09$	$8.16$	arcsec
	${\mu }_b$	$\dots$	$\dots$	$\dots$	$20.9$	$20.45$	mag arcsec-2
	$\alpha$	$\dots$	$\dots$	$\dots$	$1.50$	$2.30$
	$\beta$	$\dots$	$\dots$	$\dots$	$2.66$	$2.1$
	$\gamma$	$\dots$	$\dots$	$\dots$	$0.06$	$0.15$

Table 3. Table 3: Best-fit Imfit model for the i 𝑖 i -band image of Holm 15A. Column 1: component used in fit. Column 2: parameter. Column 3: best-fit value for parameter and 1 − σ 1 𝜎 1-\sigma confidence limits from 200 rounds of bootstrap resampling. Column 4: units. Note that for the GaussianRing3D component, we fixed the ring PA and ellipticity to both be zero, so these are not listed in the table.

Component	Parameter	Value	units
Core-Sérsic	PA	$141.9\pm 0.2$	deg
	$ϵ$	$0.187\pm 0.002$
	$n$	$0.965\pm 0.005$
	$I_b$	$20.040\pm 0.012$	mag arcsec-2
	$r_e$	$12.87\pm 0.04$	arc sec
	$r_b$	$2.57\pm 0.05$	arc sec
	$\alpha$	$12.15\pm 4.1$
	$\gamma$	$0.096\pm 0.007$
Sérsic	PA	$149.0\pm 0.1$	deg
	$ϵ$	$0.413\pm 0.003$
	$n$	$1.69\pm 0.03$
	$I_e$	$24.035\pm 0.016$	mag arcsec-2
	$r_e$	$60.67\pm 0.48$	arc sec
GaussianRing3D	PA	$52.1\pm 0.9$	deg
	inclination	$81.8\pm 1.5$	deg
	$J_0$	$1.08\pm 0.03$	counts pixel-3
	$a$	$4.37\pm 0.07$	arc sec
	$\sigma$	$1.76\pm 0.05$	arc sec
	$h_z$	$2.78\pm 0.10$	arc sec

Equations

$$\rho (r,\theta )={\rho }_{\star }(r,\theta )+M_{BH}\delta (r)+{\rho }_{DM}(r),\text{centering}\text{@add@centering}$$

$${\rho }_{DM}(r)=\frac{{\rho }_0}{{\left(1+\frac{r}{r_s}\right)}^{3-\gamma }{\left(\frac{r}{r_s}\right)}^{\gamma }},\text{centering}\text{@add@centering}$$

$$\rho_{0}=\rho_{10}\left(1+10\frac{$\mathrm{k}\mathrm{p}\mathrm{c}$}{r_{s}}\right)^{3-\gamma}\left(10\frac{$\mathrm{k}\mathrm{p}\mathrm{c}$}{r_{s}}\right)^{\gamma},\centering\@add@centering$$

$$\displaystyle\log{(M_{BH}/$\mathrm{M}_{\odot}$)}={}$$

$$+(3.29\pm 0.37)$$

$$\displaystyle\log{(M_{BH}/$\mathrm{M}_{\odot}$)}={}$$

$$+(14.19\pm 0.09).$$

$$\beta =1-\frac{{\sigma }_t^2}{{\sigma }_r^2},\text{centering}\text{@add@centering}$$

Full text

galaxies: supermassive black holes – galaxies: ETG and lenticular, cD – galaxies: evolution – galaxies: formation –stars: kinematics and dynamics – galaxies: center – clusters: individual (Abell 85)

††journal: ApJ

\savesymbol

tablenum \restoresymbolSIXtablenum

1 Introduction

Holm 15A is the brightest cluster galaxy (BCG) of Abell 85. It is a very luminous ($M_{V}=$-24.8\text{,}\mathrm{m}\mathrm{a}\mathrm{g}, Kluge et al. [2019](#bib.bib55)) early-type galaxy (ETG) with a high stellar mass of $M_{\star}\gtrsim 2\times 10^{12}\ $\mathrm{M}_{\odot}. The rotational velocity of Holm 15A is $v_{rot}\lesssim$40\text{,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s} and small compared to the velocity dispersion $\sigma\sim$350\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}. This is very common among massive ETGs (e.g Emsellem et al., 2011; Cappellari, 2016; Veale et al., 2017). Despite its high overall luminosity, Holm 15A has one of the faintest known central regions of any massive galaxy.

Figure 1 compares Holm 15A’s observed light profile with Nuker models of the centers of cored ETGs from the Lauer et al. (2007a) sample, core-Sérsic models of cored ETGs with existing dynamical models from Rusli et al. (2013a) and Thomas et al. (2016), as well as non-parametric light profiles of BCGs from Kluge et al. (2019). Evidently, at radii $r\gtrsim$30\text{,}\mathrm{k}\mathrm{p}\mathrm{c} Holm15A’s surface brightness profile is characterised by a local Sérsic index $n\gtrsim 4$, typical for massive ETGs and BCGs. Holm 15A is very bright though: only a handful of other BCGs have a higher surface brightness outside the central region ($r\gtrsim$5\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}).

It is all the more striking then how faint the center of Holm 15A is compared to ETGs from all three samples, BCG or not. Indeed, among the 88 core galaxies in the Lauer et al. (2007a) sample, the faintest center is still $\sim 0.5\$\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{/}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{2} brighter than the center of Holm 15A. Among galaxies with detailed dynamical models, the difference is even larger: $\sim 2\ $\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{/}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{2} (Rusli et al., 2013b, Thomas et al., 2016, cf. Figure 1).

Such diffuse, shallow central surface brightness regions are commonly referred to as ‘cores’ and have been observed in massive early-type galaxies (ETGs) for a long a time (e.g. Lauer, 1985; Kormendy, 1985; Faber et al., 1987). As methods for the dynamical detection of supermassive black holes (SMBHs) of ETGs have grown more sophisticated in recent years, several tight scaling relations between core properties and central black holes have been established. In particular, the most massive black holes in the local universe are expected to be found in the centers of the largest, faintest cores (e.g. Faber et al., 1997; Lauer et al., 2007a; Rusli et al., 2013a; Kormendy & Ho, 2013; Thomas et al., 2016).

The contemporary view of the formation of cores in massive ETGs is that their observed properties are best explained via so-called black hole binary ‘core scouring’. Core scouring is driven by the hardening of a SMBH binary naturally formed during dissipationless mergers between ETGs which are thought to dominate the late growth processes of massive galaxies (e.g. Khochfar & Burkert, 2003; Naab et al., 2006; Boylan-Kolchin et al., 2006; De Lucia et al., 2006; Oser et al., 2010). Gravitational slingshots eject stars on predominantly radial orbits from the center of the remnant galaxy, producing a cored central light profile (e.g. Begelman et al., 1980; Hills & Fullerton, 1980; Ebisuzaki et al., 1991; Trujillo et al., 2004; Milosavljević & Merritt, 2001; Volonteri et al., 2003; Merritt & Milosavljević, 2005; Merritt, 2006a, 2013; Rusli et al., 2013a; Rantala et al., 2018). This core-formation channel can explain the fundamental characteristics of core galaxies: (1) the observed uniform tangentially biased orbit structure in cores (Milosavljević & Merritt, 2001; Thomas et al., 2014; Rantala et al., 2018) and (2) the various core-specific scaling relations between the black hole mass, core size, size of the gravitational sphere of influence and ‘missing’ light compared to the inwards extrapolation of the steeper outer light profile (from which the core ‘breaks’; Lauer et al., 2007b; Kormendy & Bender, 2009; Kormendy & Ho, 2013; Rusli et al., 2013a; Thomas et al., 2016; Rantala et al., 2018).

From a radius of $r\sim$15\text{,}\mathrm{k}\mathrm{p}\mathrm{c} inwards down to the smallest resolved scales, the light profile of Holm 15A is almost exponential (lower panel of Fig. [1](#S1.F1)). Bonfini et al. ([2015](#bib.bib11)) and Madrid & Donzelli ([2016](#bib.bib75)) interpreted this as evidence against a large core in Holm 15A. However, as Fig. [1](#S1.F1) shows, Holm 15A fits perfectly into the homology of cored BCGs/ETGs. Hopkins et al. ([2009](#bib.bib48)) suggested that nearly exponential surface brightness profiles on kpc scales could be ubiquitous among core galaxies as a relic of merger-induced star-formation bursts in early evolutionary phases prior to the actual core formation. In their analysis, Hopkins et al. ([2009](#bib.bib48)) assumed that the sphere-of-influence of the black-hole binary is much smaller than the spatial scale relevant for these “extra-light” regions. In fact, their fits including exponential components often do not well represent the actual core region. We now know that the sizes of the cores are almost identical to the sphere-of-influence radii of the central black holes (Thomas et al., [2016](#bib.bib121)). The core of Holm 15A has a size of $3-$5\text{\,}\mathrm{k}\mathrm{p}\mathrm{c} (cf. Fig. 1, Sec. 2 and also López-Cruz et al. 2014). Hence, the expected sphere of influence is so large that it interferes with the spatial scale of potential extra-light. The only other galaxy that seems to be dominated by a nearly exponential behaviour in its entire inner region may be NGC 1600 (cf. Hopkins et al. 2009). NGC 1600 has a large sphere-of-influence radius of $1.2\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$ as well. There are many processes that influence the final inner light profile of massive galaxies, like dynamical interactions between stars and the SMBH binary, early star-formation episodes, AGN feedback etc. While these processes have been studied individually (in different levels of detail, e.g. Merritt 2006b; Hopkins et al. 2009; Teyssier et al. 2011; Martizzi et al. 2012, 2013; Choi et al. 2018; Rantala et al. 2018, 2019), we currently lack of simulations that include all these processes in a consistent manner. The black hole binary core scouring process, which is likely dominant in core formation has now been studied in great detail, including the effects of different merger histories on the stellar density profile and stellar orbits in the core (Rantala et al., 2018, 2019). Here, we use dynamical models based on new spectroscopic observations with the MUSE IFU111Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO program 099.B-0193(A). to determine the mass of the central black hole and the distribution of central stellar orbits in Holm 15A. Our goal is to shed light on possible formation scenarios for the galaxy’s extreme core.

This paper is structured as follows: Section 2 describes the new $i$-band photometry of Holm 15A obtained with the Fraunhofer Telescope at the Wendelstein Observatory, as well as additional images generated from our MUSE data. Section 3 details the MUSE spectroscopy and stellar kinematics derived from them. The dynamical models and results based on the photometry and kinematics are presented in Section 4. In Section 5 we discuss these results and their implications, in particular in view of predictions from N-body simulations. We summarize our conclusions about Holm 15A in Section6.

We use the Planck $\Lambda$CDM (Planck Collaboration et al., 2018) cosmological model, $H_{0}=67.4$, $\Omega_{M}=0.315$. The redshift of Holm 15A, z = 0.055, then corresponds to a luminosity distance of $D_{L}=$252.8\text{,}\mathrm{M}\mathrm{p}\mathrm{c} and an angular diameter distance of $D_{A}=$227.2\text{\,}\mathrm{M}\mathrm{p}\mathrm{c} ($1\arcsec=$1.10\text{,}\mathrm{k}\mathrm{p}\mathrm{c}$$).

2 PHOTOMETRY

We used two image sources for our photometric analysis of Holm 15A. The first is an $i$-band image obtained with the Fraunhofer Telescope at the Wendelstein observatory using the Wendelstein Wide Field Imager (WWFI, Kosyra et al., 2014). While a $g^{\prime}$-band image was also available, the $i$-band image had significantly better seeing (Moffat FWHM from fits to multiple stars = $0\farcs 86$ versus $1\farcs 8$ for the $g^{\prime}$-band image). The isophote analysis of this image is the basis for the 3D deprojection that we use to constrain the dynamical models (Sec.  2.1). We also used this image to analyse the core region and estimate the "missing light" in the center of Holm 15A (Sec. 2.2).

The second source is an image created from the MUSE data cube, which we used to analyse Holm 15A for the presence of dust or color gradients which could potentially affect the deprojection (Sec. 2.3, also cf. Sec. 3.2 for the spectroscopic analysis).

2.1 Wendelstein image: reduction and PSF-deconvolved light profile

Holm 15A is part of the sample of 170 local BCGs that were observed by Kluge et al. (2019) with the Wendelstein Wide Field Imager. The light profiles derived for these BCGs provide a unique photometric data base, reaching down to an unprecedented deep limiting surface brightness of $\sim 30\,\mathrm{mag/arcsec^{2}}$ in the $g^{\prime}$-band (Kluge et al., 2019, cf. Figure 1). The data cover a field of $49\arcmin{}\times{}52\arcmin{}$ (pixel size $0\farcs 2$/pixel) around Holm 15A, which corresponds to a projected area of roughly $10\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}^{2}$. The radial surface brightness profile was measured by fitting ellipses to the galaxy’s isophotes, while allowing for higher order deviations from perfect ellipses, using the code from Bender & Moellenhoff (1987). To increase the spatial resolution in the inner parts of the galaxy, the central $\sim 1\arcmin\times 1\arcmin$ of the image has been point-spread function (PSF) deconvolved using 40 iterations of the Richardson-Lucy method (Lucy, 1974). The 2D-convolution is performed on images regenerated from the previously performed isophote analysis. The radial light profile from this PSF-deconvolution is the basis of our 3D deprojection that we use to constrain the dynamical models of Holm 15A. A detailed description of the observations and data reduction can be found in Kluge et al. (2019).

2.2 Core radius and missing light of Holm 15A

The core radii of massive galaxies are typically described by either the core-break radius $r_{b}$ of a “Nuker”- (Lauer et al., 1995) or core-Sérsic profile (Graham et al., 2003; Trujillo et al., 2004), or by the ‘cusp-radius’ $r_{\gamma}$, the radius where $d\log{I}/d\log{r}=-1/2$. The cusp radius only requires that a galaxy’s light profile becomes shallow in the central parts. This is clearly the case in Holm 15A and the cusp radius is well defined: $r_{\gamma}$ = $3\farcs 7\pm 0\farcs 10$ ($4.11\pm$0.11\text{,}\mathrm{k}\mathrm{p}\mathrm{c}$$). The semi-major axis length of the corresponding isophote is $a_{\gamma}=4\farcs 1\pm 0\farcs 10$, consistent with López-Cruz et al. (2014). In contrast, the concept of a core-break radius implies – in addition to central shallowness – a distinct change of the light profile from its behaviour outside of $r_{b}$ to a different behaviour interior to $r_{b}$. As we will discuss here, the light profile of Holm 15A does not exhibit a clear and distinct change but continously flattens to the smallest observed radii.

The surface brightness distribution of Holm 15A out to $r<200\arcsec$ (or $\mu_{i}<26\mathrm{mag/arcsec^{2}}$) can be represented fairly well by the sum of two Sérsic functions, where the inner component is nearly exponential with Sérsic index $n_{1}=1.26$ and $r_{e,1}=15.81\,\mathrm{kpc}$ and the outer component follows roughly a de-Vaucouleurs profile with $n_{2}=4.21$ and $r_{e,2}=208.1\,\mathrm{kpc}$ (Kluge et al., 2019). A more complex model composed as the sum of a core-Sérsic plus a Sérsic function improves the fit in the core region slightly. The break radius of this model, $r_{b}=8\farcs 96$ (cf. model cSS in Tab. 2 of App. A.1) is roughly consistent with the radius of maximum curvature of the observed light profile. However, the Sérsic parameters of the core-Sérsic component are very different from the inner Sérsic component of the model by Kluge et al. (2019) quoted above. The “steep” Sérsic index $n_{1}=5.24$ together with the fact that $r_{e,1}<r_{b}$ undermine the intended meaning of $r_{b}$ as a "break radius" and of $n_{1}$ and $r_{e,1}$ as the local Sérsic approximation to the light outside of the core. Indeed, the corresponding Sérsic part of the model does not trace the observed light profile anywhere in the inner regions of the galaxy.

To investigate this a little further, we also tried an alternative fitting approach where we separate the determination of the core parameters from the two Sérsic components: We start by fitting the sum of two (coreless) Sérsic components to the surface brightness profile outside of the core, i.e. outside of a minimum radius $r_{\mathrm{min}}$. Then, in the second step, we repeat the fit, now including also the data inside $r_{\mathrm{min}}$ but now we only vary the core parameters in the fit, while holding the inner and outer Sérsic components $n_{1}$,$r_{e,1}$ and $n_{2}$, $r_{e,2}$, and $\mu_{e,2}$ fixed. In this way we determine the Sérsic parameters before the core parameters and force the Sérsic components to approximate the light profile outside of $r_{\mathrm{min}}$. We tried a range of different $r_{\mathrm{min}}$. Below $r_{min}<r_{\gamma}\sim 4\arcsec$ (i.e. inside the core) the inner components $n_{1}$ and $r_{e,1}$ are too much affected by the core region itself. Above $r_{\mathrm{min}}=12\arcsec$, the light profile is already so steep that we are far outside the core and the models, even after fitting the core parameters, do not provide good fits anymore. For $r_{\mathrm{min}}=4\arcsec-12\arcsec$, these two-step fits represent the data very well. Moreover, in all the two-step fits we found $r_{b}<r_{e}$, and $r_{b}\sim r_{\gamma}$, as expected (cf. models cSS($r_{\mathrm{min}}=4$) and cSS($r_{\mathrm{min}}=12$) in Tab. 2). The Sérsic components approach the Sérsic+Sérsic model of Kluge et al. (2019) in the limit of small $r_{\mathrm{min}}$. However, the fits did not converge to a stable set of parameters. We found both the Sérsic index $n_{1}$ of the inner component and $r_{b}$ to systematically increase with $r_{\mathrm{min}}$ (cf. Tab. 2 and Fig. 2).

All these results led us to conclude that the galaxy does not exhibit a clear break radius inside of which the light profile follows a power law and outside of which it can be characterised by a single local Sérsic index $n$ over a range that is more extended than a few arcseconds. Fitting the inner parts of the 1D light profile of Holm 15A with a Nuker profile confirmed this finding. Again, we could not derive a stable break radius and $r_{b}$ turned out to be a monotonic function of the maximum radius out to which we extended the fit (we tried $r_{\mathrm{max}}=10-70\,\arcsec$; cf. Tab. 2).

Finally, we also performed a 2D- multi-component fit to the entire $i$-band Wendelstein image of Holm 15A using IMFIT (Erwin 2015; see Appendix A.2). This yielded a stable set of core parameters. However, in the 2D-analysis, allowing for a broken inner profile with a power-law core did not improve the fit significantly over a central, pure Sérsic component with $n\sim 1$ and $r_{b}=2\farcs 57$.

Holm 15A evidently continues the homology of cores observed in less extreme ellipticals in the sense of having a faint center with a shallow surface brightness profile (cf. Fig. 1). But, as our attempts of identifying a clear break radius have shown, the core region in Holm 15A is not as sharply separated from the outer parts of the galaxy as it is in other core galaxies with a more prominent break in the light profile. Because of this, even though both $r_{b}$ and $r_{\gamma}$ have been shown to follow tight scaling relations with $M_{BH}$ in other core galaxies (e.g. Lauer et al., 2007b; Thomas et al., 2016), we will only consider the cusp radius of Holm 15A in the rest of the paper.

The shallowness of the inner light profile still allows the estimation of the amount of "missing light". From the above described models cSS($r_{\mathrm{min}}=4$) and cSS($r_{\mathrm{min}}=12$) (see Tab. 2) we find $L_{i,def}=(2.75\pm 2.22)\times 10^{10}L_{i,\odot}$, which we will later use in Section 5.1 to estimate the mass of stars ejected from the center via core scouring. The estimated missing light is illustrated in Fig. 2.

2.3 MUSE images: no evidence for dust or color gradients

To investigate whether dust extinction might distort the isophotes, and to check for color gradients indicative of a change in the stellar populations, we also generated images from the MUSE data cube. This has two advantages. First, the MUSE observations have (slightly) better seeing than the Wendelstein $i$-band image: in the “red” image (see below for definition), we measured FWHM = 0$\farcs$72 from the two point sources in the image. Second, when collapsing the data cube we can choose wavelength ranges that explicitly exclude emission, which is important because we do detect regions of line emission within Holm 15A (see below).

We use the spectral region 7300–8500 Å to create a largely emission-line-free “red” image and the spectral region 4750–5500 Å for its “blue” counterpart. The ratio of the blue and red MUSE images is shown in the right-hand panel of Figure 3 and it shows no evidence for either dust lanes or significant color gradients.

2.4 3D deprojection

In order to constrain the distribution of stars in our dynamical model of Holm 15A (see Section 4), we create a 3D deprojection of the luminosity density from our deconvolved 2D Wendelstein image. The algorithm that we use to achieve this enables us to find a 3D non-parametric axisymmetric luminosity density distribution $\nu(\bm{r})$ consistent with the 2D input surface brightness distribution and an assumed inclination angle $i$. As can be seen in Figure 2, Holm 15A is for the most part relatively round, but flattens significantly to an ellipticity $\varepsilon\sim 0.4$ at radii $\gtrsim 100\arcsec$. In the axisymmetric case, this limits possible viewing angles to be close to edge on, which is why we assume $i=$$$. The algorithm utilizes a penalized log-likelihood function and is detailed in Magorrian (1999). As Figure 2 shows, the resulting axisymmetric luminosity density distribution reproduces the relevant observed photometric features almost perfectly.

3 MUSE spectroscopy: stellar kinematics of Holm 15A

3.1 MUSE observations and data reduction

We obtained wide-field spectroscopic data of Holm 15A from the Multi-Unit Spectroscopic Explorer (MUSE) at the Very Large Telescope at Paranal on 2017 November 16 and 2018 August 10. At z = 0.055 MUSE covers several important absorption features such as H$\beta$, the Mgb region, NaI, several Fe absorption features and the Ca II triplet.

Our observations were carried out over the course of two nights and consist of three observational blocks of two dithered $1200\text{\,}\mathrm{s}$ exposures of Holm 15A plus one $300\text{\,}\mathrm{s}$-long exposure of the sky, inbetween each. All observations, including the sky-field offset, cover an approximately 1\text{\,}\mathrm{\SIUnitSymbolArcminute}$\times$1\text{\,}\mathrm{\SIUnitSymbolArcminute} FOV composed of 24 combined integral field units (IFUs).

We performed the data reduction using version 2.8.5 of the standard Esoreflex MUSE pipeline supplied by ESO (Freudling et al., 2013). The pipeline runs several recipes on both exposures such as flat-field and wavelength calibrations and returns a combined data cube, covering the optical domain from about $4800\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ to $9400\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ with a spectral resolution of $1.25\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. We sampled the cube in spaxels of $0\farcs 4\times 0\farcs 4$, which at the redshift of the galaxy (z = 0.055) corresponds to approximately 400\text{\,}\mathrm{p}\mathrm{c}$\times$400\text{\,}\mathrm{p}\mathrm{c} per pixel. As previously mentioned, we measure a PSF with FWHM = $0\farcs 71$ for the MUSE image.

Sky emissions were removed separately from all galaxy exposures using the sky-field from offset sky-exposures, taking into account the instrumental line spread function for each IFU.

3.2 Treament of spectra and derivation of (parametric) stellar

kinematics

For our study of Holm 15A, we initially used the MUSE absorption spectra to derive spatially resolved, 2D stellar kinematics parameterized by the rotational velocity $v_{rot}$, velocity dispersion $\sigma$ and higher-order Gauss-Hermite coefficients $h_{3}$ and $h_{4}$ of the line-of-sight velocity distribution (LOSVD). For the dynamical modelling, we use non-parametric LOSVDs that were derived following a set of equivalent steps (see Sec. 3.3).

To achieve a balance between a precise measure of the kinematics in the core and an overall high spatial resolution we aim for a target S/N of at least $\sim 50$ per pixel in each spectrum. To achieve this, we spatially bin the data cube using the Voronoi tessellation method of Cappellari & Copin (2003). Pixels belonging to foreground sources such as galaxies or AGN are removed from the data before binning.

At the center of the galaxy ($r\leq 5\,\mathrm{kpc}$) the spatial resolution of the Voronoi bins turns out to be $0.4^{\prime\prime}$ - $0.8^{\prime\prime}$ (roughly $400-$800\text{,}\mathrm{p}\mathrm{c}) for a S/N $\sim 50$. We here define the radius of the gravitational sphere of influence (SOI) of the black hole as the radius where the enclosed mass $M(\leq r_{SOI})\equiv M_{BH}$. By integrating the deprojected 3D luminosity density and assuming a range of plausible stellar mass-to-light ratios, between $\Upsilon_{*}=4$ and $6$, we estimated the enclosed mass of the galaxy. For the lowest expected black hole mass for a galaxy of this mass and velocity dispersion, $M_{BH}\sim 3\times 10^{9}$\mathrm{M}_{\odot} (using the mean expected values from the $M_{BH}-\sigma,M_{Bu}$ scaling relations for ETGs from McConnell & Ma, 2013; Saglia et al., 2016)), the enclosed stellar mass equals $M_{BH}$ at $r_{SOI}\sim 1\farcs 6$. Since our PSF and spatial binning resolution are both on the order $0\farcs 8$ we ensure that we can resolve the expected sphere of influence (SOI) with a diameter of $2\times 1\farcs 6=3\farcs 2$ by a factor $\geq 4$. However, the extreme core properties of Holm 15A actually point to a SMBH with $M_{BH}\sim 10^{11}M_{\odot}$ (based on $M_{BH}-r_{\gamma}$ scaling relations from Lauer et al. (2007c) and Thomas et al. (2016)), whose SOI radius would be roughly $r_{\mathrm{SOI}}\sim 4-5\,\arcsec$ – a factor $>10$ above our resolution limit. If the dark matter halo is included in the modeling, this resolution is sufficient for a robust black hole mass determination (Rusli et al., 2013b).

In total, we obtain 421 spatial bins, of which 145 bins are located inside the central $5^{\prime\prime}$. For the purpose of our subsequent dynamical modeling of the galaxy we divided the spatial bins of our MUSE FOV into four quadrants, q1-4 in such a way that quadrant membership is determined by which side of the major and minor axes the center of each bin is located on

Parametric LOSVDs for each bin were obtained by fitting the stellar absorption lines of the galaxy with Penalized Pixel-Fitting (pPXF, Cappellari, 2017) implemented in Python 2.7. PPXF convolves a weighted sum of template stellar spectra, in this case the MILES library (Sánchez-Blázquez et al., 2006) with a Gauss-Hermite LOSVD in order to fit the absorption features. Optionally, emission-line features of ionized gas are fit simultaneously, with a separate set of templates and LOSVDs. Figure 4 shows an example of a (parametric) kinematic fit to the spectral features of Holm 15A with pPXF for a bin located roughly $0\farcs 5$ from the center of the galaxy (best fit to stellar component: $v_{rot}=-1.59\pm$8.04\text{,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s} relative to the systemic velocity of the galaxy, $\sigma=342\pm$9.71\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}, $h_{3}=0.025\pm 0.015$, $h_{4}=0.062\pm 0.018$).

Several bins within the central $5\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$ of the galaxy – primarily in the southeastern regions – region contain emission lines from ionized gas, most notably H$\alpha$, H$\beta$ , [OIII] $5007\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, [NI] $5199\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ and [NII] $6583\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ (cf. Figure 4), which we fitted with the emission line fitting routine of pPXF, though we do not consider their kinematics in this study. Figure 5 shows the measured emission line flux for H$\alpha$, H$\beta$, [OIII] and [NII]. The average flux ratios $\log(\mathrm{[OIII]/H\beta})=0.09\pm 0.26$ and $\log(\mathrm{[NII]/H\alpha})=0.48\pm 0.12$ of emission lines with S/N > 3 are associated with LINER-type emission (Kauffmann et al., 2003), which is quite typical for cool-core clusters. Of the $\sim 100$ brightest X-ray clusters, Abell 85’s cool core has the 14th strongest cooling flow (Chen et al., 2007). The spatial extent of this LINER-type emission ($\sim 4$–$5\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$) suggests it could be related to ionized cooling-flow filaments (e.g Ferland et al., 2008, 2009; Ogrean et al., 2010). This was already previously noted by McDonald et al. (2010), who found that it coincided with a similarly extended region of X-ray emission associated with cooling flows.

By contaminating some absorption features such as H$\beta$, the gas emission increase the uncertainties of the kinematic fits in some bins. As we will show in section 4, this contamination of mostly central spectra slightly increases the uncertainty of $M_{BH}$, but has little impact on the global stellar mass-to-light ratio $\Upsilon_{*}$ and the shape of the dark matter halo.

At redshift z = 0.055, the strong oxygen $5577\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ sky emission line lies on top of the $5270\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ Fe-feature. Because this line is difficult to remove, the Esoreflex sky subtraction left strong residuals in this region, effectively rendering it unusable for fitting. We noted a few additional systematic residuals which may be related to sky subtraction or telluric correction issues as well. In order to minimse possible systematics in the LOSVDs, we defined a single mask that we used for all spectra throughout the entire galaxy. We consistently mask all wavelength regions that are possibly affected by any systematic issues.

We performed our kinematic fits over the spectral interval between 5010 and $7050\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. Including spectral regions bluer than $5010\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ resulted in lower-quality fits and a constant bias in $h_{3}$, indicative of template mismatch. Spectral regions redder than $7050\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ were badly affected by sky lines and were therefore omitted. In particular, we could not derive meaningful kinematics in the [Ca II] triplet region.

We also used a 6th order multiplicative polynomial and an additive constant in the fit. The former allows for the correction of errors in the flux calibration, while the latter is typically used to correct over- or underestimations of the continuum during sky correction. We also made use of the sigma clipping and bias-factor options. The value of the bias factor – 0.2 in our case – was determined from testing pPXF on Monte Carlo simulations of model spectra.

A subset of stellar template spectra for the fit was selected as follows: We fitted a mean spectrum of all bins of the galaxy with the full set of 985 MILES library templates. All binned spectra were corrected for the systematic velocity of the galaxy, as well as their respective rotational velocities. All spectra were normalized to one before averaging. We set both the third-order Gauss-Hermite coefficient $h_{3}$ and the additive constant to zero in order to avoid template mismatch (which can result in biases in these parameters). With these restrictions pPXF assigned non-zero weights exclusively to a set of 16 templates with a wide variety of luminosity classes but limited to spectral types G, K and M, in good agreement with the uniformly red color of the galaxy (ec. 2). We used this subset of stars from the MILES library as templates for fitting the galaxy’s absorption features in all Voronoi bins.

The parameterized kinematics in the interval between 5010 and $7050\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ over the MUSE FOV are shown in Figure 6. As can be seen in the figure, we measure a weak rotation signal of less than $40\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$, which is only faintly reciprocated in $h_{3}$ – the rotation is likely too weak for an anti-correlated signal in this parameter to be detectable. The velocity dispersion $\sigma$ peaks in the central regions ($r<$2\text{,}\mathrm{k}\mathrm{p}\mathrm{c}) at $\sim$350\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}, stays somewhat constant at $\sim$ $330\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$ throughout most of the FOV and finally starts to rise again at the edges of the MUSE FOV up to $\gtrsim$370\text{,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$. Our measured velocity dispersions are similar to those of Fogarty et al. (2014). Our $h_{4}$ kinematic profile starts out at $\sim 0.07$ within $2\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$ and rises to $\gtrsim 0.1$ towards the edges of the FOV. In Appendix B.1 we compare the kinematics of Holm 15A to those of massive ETGs from the MASSIVE survey. The corresponding statistical uncertainties are shown in Figure 7. Uncertainties were determined from Monte Carlo simulations on model spectra of the galaxy, i.e. re-fitting best-fit spectral models with 100 different noise realizations, the noise being drawn from a Gaussian distribution with a dispersion corresponding to the local S/N, which is measured directly from each spectrum. We note that the distribution of uncertainties is spatially asymmetric between central bins across quadrants – central kinematics in q3 have overall larger uncertainties than those in the other quadrants. This is in agreement with the distribution of emission-line flux between quadrants (cf. Figure 5), i.e. q3 seems to be affected worse by uncertainties introduced by gas contamination of absorption features. However, as we will show in Section 4, including q3 in our dynamical modeling did not produce any larger systematic offset in our best fit parameters relative to the other quadrants.

3.3 Non-Parametric LOSVDs

In our dynamical modeling of Holm 15A we set out to achieve a precise mass measurement of the galaxy, which makes the parametric representation of the stellar kinematics in Figure 6 problematic: Large values of $\sigma$ and $h_{4}>0$ over the entire FOV result in the escape velocity of the galaxy, $v_{esc}$ being practically infinite everywhere. Since $v_{esc}$ depends directly on the gravitational potential we try to measure it as accurately as possible.

To obtain LOSVDs with more realistic $v_{esc}$, we use our own kinematic extraction code (Thomas et al. in prep.) which operates in a similar way as pPXF but minimizes the $\chi^{2}$ over all spectral pixels by utilizing a Levenberg-Marquardt algorithm to fit a template broadened with a non-parametric LOSVD to the absorption features of a galaxy.

We use the same setup of template stars, additive and multiplicative polynomials as described above. Emission lines are masked for each spectrum individually, according to their respective widths (spectral regions within $4\times\sigma_{gas}$ are masked for each emission line) and positions as determined with the pPXF emission line fit. The non-parametric LOSVDs mainly differ from the parametric ones in the high-velocity tails, as demonstrated for an example bin of Holm 15A in Figure 8. While the width of the LOSVD ($\sigma=338\pm$9.57\text{,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s} with our own code and $\sigma=328\pm$10.7\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s} with pPXF), as well as its global shape, are similar for both methods, the non-parametric LOSVDs provide a more realistic sampling of the LOSVD and noise at large projected velocities. Therefore, for our dynamical study of Holm 15A, we use the non-parametric LOSVDs. Radial profiles comparing both parametric and non-parametric kinematics for all bins in our study are presented in Appendix B.2.

4 Schwarzschild Dynamical Modeling of Holm 15A

4.1 Dynamical models

We dynamically modeled Holm 15A under the assumption of axisymmetry. The lack of unambiguous, obvious isophotal distortions (see Section 2) and the overall symmetry of the observed kinematic profiles (see Section 3.2) imply that Holm 15A is generally consistent with an axially symmetric stellar distribution.

The dynamical models in this study were constructed using an updated version of our axisymmetric Schwarzschild orbital superposition code. We will here only briefly summarize the key features of our implementation and refer to previous publications for more in-depth descriptions (Richstone & Tremaine, 1988; Gebhardt et al., 2003; Thomas et al., 2004; Siopis et al., 2009).

Schwarzschild dynamical modeling is based on the calculation of stellar orbital distributions in a fixed gravitational potential as a solution to the collisionless Boltzmann equation (Schwarzschild, 1979). Any orbit can be fully described by three integrals of motion: The classical integrals $E$ and $L_{z}$ (in the axisymmetric case) plus a non-classical integral $I_{3}$ (in most astrophysically relevant cases). Sampling values of this set of integrals of motion ($E$, $L_{z}$, $I_{3}$) allows us to create an orbit library in a given gravitational potential $\Phi$ whose distribution function $f(\bm{r},\bm{v})$ satisfies the collisionless Boltzmann equation.

In order to determine $\Phi$, we assume that the density distribution of Holm 15A can be described by

[TABLE]

which we insert into Poisson’s equation. $\rho_{\star}$ is linked to the three dimensional deprojection $\nu(r,\theta)$ of the observed i-band surface brightness (cf. Section 2) via the stellar (i-band) mass-to-light ratio, $\rho_{\star}(r,\theta)=\Upsilon_{\star}\cdot\nu(r,\theta)$, assuming a spatially constant stellar $\Upsilon_{\star}$. In addition to the mass of the central black hole $M_{BH}$, the model admits the inclusion of a dark matter (DM) halo $\rho_{DM}(r)$. Here, we chose a generalised NFW-halo derived from cosmological N-body simulations (Navarro et al., 1996; Zhao, 1996):

[TABLE]

with

[TABLE]

where $\rho_{10}$ is the DM density at $10\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$, $r_{s}$ the scale radius of the halo and $\gamma$ the inner slope of the DM density profile.

For a given $\Phi$, we sample thousands of representative initial orbital conditions, implicitly varying all the integrals of motion $E$, $L_{z}$ and $I_{3}$, and including individual orbital phase-space volumes (Thomas et al., 2004). For Holm 15A, we stored LOSVDs in 29 velocity bins adapted to the velocity dispersion of the galaxy, with one LOSVD associated with each of the 421 spatial bins of our FOV, meaning our models fitted roughly a total of 3000 velocity bins per quadrant.

We use the NOMAD optimization software (Audet & Dennis, Jr., 2006; Le Digabel, 2011; Audet & Hare, 2017) to find the set of mass parameters $M_{BH}$, $\Upsilon_{\star}$, $\rho_{10}$, $r_{s}$ and $\gamma$ that yields the best fit to the observed kinematics.

4.2 Results

The most important result from our dynamical modeling is the detection of a SMBH with $M_{BH}=(4.0\pm 0.80)\times 10^{10}$\mathrm{M}_{\odot} in Holm 15A. The associated SOI of this SMBH is $r_{SOI}=3.8\pm$0.37\text{\,}\mathrm{k}\mathrm{p}\mathrm{c} ($3\farcs 5\pm 0\farcs 34$). Even though the galaxy is more than 200 Mpc away, we spatially resolve the SOI by a factor of 10. In fact, $\sim 100$ out of our 421 LOSVDs sample the SOI of the galaxy. The modeling results for the black hole, stellar mass-to-light ratio and DM halo parameters are summarised in Table  1. $\Delta\chi^{2}$ curves for $M_{BH}$, $\Upsilon_{\star}$ and $\rho_{10}$ from all four quadrants are shown in Figure 9. The figure shows that none of the four quadrants stands out and yields a significantly different result than the others. While the black hole mass in q3 (where the gas emission in the spectra is most prominent) is slightly larger than in the other quadrants, this offset is not significant. By computing the dynamical quantities separately for each quadrant and estimating the uncertainties from these four nearly independent measurements, we implicitly include any residual systematics (like, e.g., from the gas emission) in our error budget. Fits to the kinematics of one quadrant of Holm 15A parameterized by $v_{rot}$, $\sigma$, $h3$ and $h4$ of our best-fit model are shown in Figure 10. They show that our best-fit model can successfully reproduce the observed kinematics of the galaxy. For the non-parametric kinematics our best-fit model reaches a reduced $\chi^{2}$ of $0.8-0.9$ for each quadrant.

We had previously also acquired spectroscopy of Holm 15A from the McDonald Observatory using the low-resolution mode ($\sigma\sim$25\text{,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$) of the integral field unit spectrograph VIRUS-W (Fabricius et al., 2012). Stellar kinematics for these independent data were derived by applying the Fourier Correlation method (FCQ) by Bender (1990) in the wavelength interval between 4500 and 6250 Å, using a sparser spatial sampling (Figure 10, blue) and circular spatial binning. This entirely independent measurement of the stellar kinematics in Holm 15A is consistent with the MUSE kinematics. We note that on average values of $h_{4}$ and $\sigma$ appear to be slightly lower for FCQ (likely due to a different smoothing-method). Therefore, as a consistency check we ran a second set of dynamical models using only the VIRUS-W kinematics and found the same results within the errors. Because the MUSE data have better spatial resolution and higher signal-to-noise ratio we will only discuss the results derived from the MUSE data in the remainder of this paper.

Finally, an example comparison between an observed and modelled LOSVD and a discussion of the importance of the LOSVD wings can be found in Appendix C.

5 Discussion

With $M_{BH}=(4.0\pm 0.80)\times 10^{10}$\mathrm{M}{\odot}, the SMBH at the center of Holm 15A is the most massive dynamically determined black hole so far. It is a factor of two larger than the SMBHs in NGC 4889 (McConnell et al., [2012a](#bib.bib81)), with $M_{BH}=(2.1\pm 0.99)\times 10^{10}$\mathrm{M}_{\odot} and NGC 1600 with $M_{BH}=(1.7\pm 0.15)\times 10^{10}$\mathrm{M}{\odot} (Thomas et al., [2016](#bib.bib121)). Quasar luminosities at higher redshifts and current determinations of local SMBH scaling relations give an expected black hole cumulative space density ranging from half a dozen up to a few hundred SMBHs with $M_{BH}\gtrsim 10^{10}$\mathrm{M}_{\odot} out to $z\leq 0.055$ (e.g. Lauer et al., 2007b; Rusli et al., 2013b). Hence, circumstances for the formation of a 40-billion-solar-mass SMBH are probably rare, but the central structure of the Coma cluster serves as an example that they do exist. As stated above, NGC4889, one of the two central galaxies of Coma, contains a SMBH of $M_{BH}=2.1\times 10^{10}$\mathrm{M}_{\odot}$$. The other galaxy, NGC4874, has a very extended classical shallow-power-law surface-brightness core with a size of $r_{b}=1.7\,\mathrm{kpc}$ (Lauer et al., 2007a). This suggests a SMBH with a mass of $M_{\mathrm{BH}}\sim 2\times 10^{10}\,M_{\odot}$ (using the core scaling relations of Thomas et al. 2016). Both galaxies are in interaction and will eventually merge (e.g. Arnaboldi et al., 2006; Gerhard et al., 2007). This will produce a BCG at the center of the Coma cluster which will very likely have a SMBH in the same mass range as Holm 15A has now.

In the following sections we will discuss the observational and theoretical evidence for the merger origin of Holm 15A, as well as attempt to unravel some specific details of the merger history.

5.1 SMBH-scaling relations: Evidence for dissipationless merging

The SMBH of Holm 15A is not only the most massive one to date, it is also four to nine times larger than expected given the galaxy’s stellar mass $M_{Bu}=(2.5\pm 0.64)\times 10^{12}$\mathrm{M}_{\odot}$$ and the galaxy’s stellar velocity dispersion $\sigma=(346\pm 12.5)\ km/s$ (see Figure 11a).

It has been previously noted that the $M_{BH}-\sigma$ relation may shallow out at the high-mass end due to dry merging becoming the dominant growth process at the high-mass end. Since dry (major) mergers grow $\sigma$ only slowly (e.g. Lauer et al., 2007b; Naab et al., 2009) but simply sum over the central SMBH masses of the merging galaxies, such mergers will move galaxies towards “overmassive” $M_{BH}$ at a given $\sigma$ (e.g. Lauer et al., 2007b; Kormendy & Bender, 2013). Correspondingly, massive core galaxies follow a $M_{BH}-\sigma$ relation that is steeper and slightly offset (towards larger values of $M_{BH}$) compared to less massive, cuspy galaxies (cf. Saglia et al. 2016 and McConnell & Ma 2013). Despite the fact that we here already consider the $M_{\mathrm{BH}}-\sigma$ relation of core galaxies, Holm 15A is still almost an order of magnitude offset in $M_{BH}$ (see Figure 11a). This might be indicative of an especially extensive dry merging period.

One could expect the $M_{BH}-M_{Bu}$ relation to be tighter at the high-mass end, since the ratio $M_{\mathrm{BH}}/M_{\mathrm{Bu}}$ ratio is conserved in dry mergers. Holm 15A, however, is also a strong outlier from this relation ($M_{BH}$ is roughly 4 times larger than expected from $M_{Bu}$, see Figure 11b). The ratio between $M_{BH}$ and $M_{Bu}$ is typically $\lesssim 0.5\%$ for cored ETG and typically $\lesssim 1\%$ when considering all ETGs below a stellar mass of $<10^{13}M_{\odot}$, irrespective of central morphology (Kormendy & Ho, 2013). Holm 15A, however, hosts a black hole that contains close to $2\%$ of the total stellar mass of the galaxy. A similar high ratio as been found in NGC 1600 (Thomas et al., 2016). This might suggest that the progenitor galaxies of Holm 15A were different from typical massive ETGs at $z\sim 0$. Studies of the evolution of $M_{\mathrm{BH}}/M_{\mathrm{Bu}}$ since $z\sim 3$ in active galaxies suggest that the ratio scales like $(1+z)^{0.7-1.4}$ (e.g. Decarli et al., 2010; Merloni et al., 2010; Bennert et al., 2011). Depending on which $M_{BH}-M_{Bu}$ relation is used (all central morphologies or cores-only) we can estimate that Holm15A’s progenitors might have formed early, at $z\gtrsim 1$ or $2$.

In Fig. 11b we only consider scaling relations based on dynamcial bulge masses to avoid systematics related to assumptions about the initial stellar mass function (IMF). We will touch on this again in Sec. 5.6.

Kluge et al. (2019) showed that BCGs and ETGs in general follow different scaling relations between total luminosity, size and effective surface brightness. This would also translate into different SMBH scaling relations. Bogdán et al. (2018) suggested that BCGs follow steeper $M_{BH}-\sigma$ and $M_{BH}-M_{Bu}$ relations (cf. Fig. 11a,b). Holm 15A is closer to these BCG-centric scaling relations. In fact, it happens to fall onto the corresponding $M_{BH}-M_{Bu}$ relations and is offset from the corresponding $M_{BH}-\sigma$ relations by about a factor of two. This could indicate that the galaxy formed from a dissipationless, (roughly) equal-mass BCG-merger, though the scatter in the relations is large.

We note that the total stellar mass of Holm 15A is estimated based on the assumption that the mass-to-light ratio is constant out to a region that is almost 10 times larger than the field of view of our kinematic observations. Therefore, in Fig. 11c, we also compare Holm 15A’s K-band luminosity $L_{K}$ to the $M_{BH}-L_{K}$ relation of Kormendy & Ho (2013). $L_{K}$ was measured from an image that extends out to $\sim$250\text{,}\mathrm{k}\mathrm{p}\mathrm{c}$$ and that was obtained with the three-channel imager at the Wendelstein 2-m Telescope (3KK Lang-Bardl et al., 2010, 2016). Holm 15A follows the $M_{BH}-L_{K}$ correlation better than the $M_{BH}-M_{Bu}$ relation.

5.2 Scaling relations of core properties: similarity with other core galaxies

Dissipationless mergers between ETGs involve binary black hole core-scouring and, hence, result in depleted, low-surface-brightness cores. As already mentioned above, core galaxies follow specific scaling relations between the core size, missing light and black-hole mass (Lauer et al., 2007a; Rusli et al., 2013a) and the radius of the sphere-of-influence (Thomas et al., 2016). Simulations have shown that these relations can be explained by the black-hole binary model (Rantala et al., 2018). In Figure 11d we show the central stellar mass deficits from Rusli et al. (2013a) together with Holm 15A. The mass deficit in Holm 15A is $M_{\star,def}=(1.24\pm 1.00)\times 10^{11}M_{\odot}$, based on the dynamical stellar mass-to-light ratio and $L_{i,def}$ derived in Sec. 2.2 . This roughly corresponds to $0.5-5.5$ times the black hole mass, similar to the mass deficits in many other core galaxies.

In Figure 11e we compare the core size of Holm 15A to other galaxies. As described in Sec. 2.2 we use the cusp radius $r_{\gamma}$ here. Compared to the galaxies of Rusli et al. (2013b, a); Thomas et al. (2016) the core in Holm 15A is roughly a factor $2.5$ larger than expected for the mass of its black hole.

Such an offset could be explained, for example, if Holm 15A experienced an early phase of rapid evolution with an enhanced merger rate. It could well be then that not only a binary black hole was involved in the formation of its core, but possibly a more complicated system of multiple black holes. Theory suggests that core scouring efficiency is significantly enhanced by multiple black holes and that cores grow much larger (Kulkarni & Loeb, 2012). We will revisit this issue in Sec. 5.4. In Figure 11f we compare $r_{\gamma}$ with the radius of the sphere of influence $r_{SOI}$. Despite being offset on the $M_{BH}-r_{\gamma}$ relation the cusp radius is consistent with the correlations between core-size measurements and $r_{SOI}$ in other core galaxies

5.3 A new correlation between black-hole mass and core surface brightness

Cores in massive ETGs obey a strong homology in that the central surface brightness correlates inversely with the size of the core (Faber et al., 1997; Lauer et al., 2007b) – This, together with the scaling between $M_{BH}$ and core size, implies a potential scaling between $M_{BH}$ and the central surface brightness $\mu_{0}$ in cores. An equivalent argument can be made for a correlation between $M_{BH}$ and central stellar surface mass density $\Sigma_{0}$. We show these correlation in Figure 12 for the galaxy sample of Rusli et al. (2013b), NGC1600 (Thomas et al., 2016) and Holm 15A. We used the uncertainties for the stellar mass-to-light ratios and black hole masses listed in Rusli et al. (2013a, b) and (Thomas et al., 2016) and assumed rather conservative uncertainties of $0.1\text{\,}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{/}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{2}$ for the light profiles. Our best-fit linear relations were determined following the approach to linear regression from Kelly (2007) (using the Python package linmix by Meyers 2015) with errors in both $M_{BH}$ and $\mu_{V,0},\Sigma_{0}$:

[TABLE]

[TABLE]

The $M_{BH}-\mu_{V,0}$ relation has an intrinsic scatter $\epsilon=0.32\pm 0.07$. Similarly, the $M_{BH}-\Sigma_{0}$ relation has an intrinsic scatter of $0.30\pm 0.07$. Values of $\Sigma_{0}$ were calculated from the surface brightness at the spatial resolution limit for each galaxy and their corresponding dynamical stellar mass-to-light ratios (Rusli et al., 2013a, b; Thomas et al., 2016). Values for both $\mu_{V,0}$ and $\Sigma_{0}$ were determined using the observed light profiles of each core galaxy. Holm 15A has the lowest central stellar surface brightness/mass, $\mu_{V,0}$ = $19.9\pm$0.13\text{,}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{/}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}^{2}, $\Sigma_{0}=(3.0\pm 0.40)\times 10^{3}$\mathrm{M}_{\odot}\mathrm{/}\mathrm{p}\mathrm{c}^{2} of all core galaxies with dynamical black hole mass measurements (cf. Figure 12). Nonetheless, Holm 15A is fully consistent with the homology established by other core galaxies. 222The listed relations were determined including Holm 15A, but the relations change only marginally and within the listed uncertainties when we exclude the galaxy. All of the above evidence points to the fact that the core in Holm 15A was formed by the same physical process as cores in other massive ETGs, i.e. by a black-hole binary.

5.4 N-body merger simulations: evidence for a merger between two core galaxies

We will now discuss what the specific photometric and orbit-dynamical properties of Holm 15A may tell us about its merger history.

In Figure 13 we compare the light profile of Holm 15A with the N-body merger simulations of Rantala et al. (2018, 2019). These simulations study the outcome of a dissipationless merger between two early-type progenitor galaxies, both with central black holes. The simulations follow the dynamical interaction between the black hole binary that temporarily forms at the center of the remnant galaxy and the surrounding stars with high accuracy. The figure demonstrates that mergers between cuspy progenitors (i.e. mergers between originally coreless progenitor galaxies) lead to slightly different light profiles than do mergers between galaxies that already had cores. The light profile of Holm 15A, in fact, looks very similar to the 2nd type of merger, i.e. between two already cored galaxies333At roughly $8\times r_{\gamma}\sim$40\text{,}\mathrm{k}\mathrm{p}\mathrm{c} (for Holm 15A) the surface brightness of the rescaled core-core remnant drops faster than the that Holm 15A. This could be due to the fact that the merger simulations do not include an extended cD halo. Photometric studies of Holm 15A (e.g. Kluge et al., [2019](#bib.bib55); Donzelli et al., [2011](#bib.bib25)) suggest an extended stellar envelope starting at $r\gtrsim$35\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}. At radii ¡ $8\times r_{\gamma}$ the core-core remnant is remarkably similar to Holm 15A. (Figure 13).

The evidence in favor of a core-core merger from the light profile is consistent with the evidence from the orbit distribution that we find in Holm 15A. Figure 14 shows the radial profile of the anisotropy parameter

[TABLE]

where $\sigma_{r}$ is the radial and $\sigma_{t}=\sqrt{(\sigma_{\theta}^{2}+\sigma_{\phi}^{2})/2}$ is the tangential velocity dispersion, computed from the dispersions $\sigma_{\theta}$ and $\sigma_{\phi}$ in the two angular directions. The figure also includes the results from the numerical N-body simulations. It is known that core scouring results in an orbital distribution that is biased increasingly towards tangential orbits ($\beta<0$) inside the SOI of the black hole as $r\rightarrow 0$ and increasingly towards radial orbits ($\beta>0$) outside of it, towards larger radii (e.g. Quinlan & Hernquist, 1997; Milosavljević & Merritt, 2001; Rantala et al., 2018). Tangential anisotropy around SMBHs has been observed in systems of various masses and morphologies (e.g. Verolme et al., 2002; Gebhardt et al., 2003; Shapiro et al., 2006; Houghton et al., 2006; Gebhardt & Thomas, 2009; Gültekin et al., 2009; Krajnović et al., 2009; Siopis et al., 2009; Shen & Gebhardt, 2010; van den Bosch & de Zeeuw, 2010; Schulze & Gebhardt, 2011; Gebhardt et al., 2011; McConnell et al., 2012b; Walsh et al., 2015; Feldmeier-Krause et al., 2017; Thomas et al., 2016). In core galaxies, specifically, the measured anisotropy is extremely homogeneous and intimately linked to the core region and follows very closely the prediction of N-body merger simulations (Thomas et al., 2014).

In Holm 15A we see the same behaviour: a change from outer radial anisotropy to inner tangential motions roughly at the sphere of influence radius (which is similar to the core size, see Figure 11f). The evidence for this comes from the wings of the observed LOSVDs (cf. App. C). However, the central anisotropy in Holm 15A is milder than observed in other core galaxies, which follow the “cuspy-cuspy” line in Figure 14 (Rantala et al., 2018). This difference is actually expected if the direct progenitors of Holm 15A were not cuspy power-law ellipticals but galaxies that already had cores. In the latter case, the anisotropy in the center is predicted to be very similar to the observed orbital structure of Holm 15A (Rantala et al., 2019)444In the N-body simulations, the final anisotropy profile of an equal-mass core-core merger is very similar to that of the final orbit distribution after a sequence of minor mergers (Rantala et al., 2019). However, the light profile of Holm 15A is more similar to the core-core merger than to the remnant after repeated minor mergers. Further simulations covering a wider range of initial conditions are needed to confirm the connection between anisotropy, profile shape and merger history..

Since cores grow with each merger generation, a core-core merger scenario would plausibly explain the fact that the central region of Holm 15A is fainter than the centers of $\gtrsim 97\%$ of the 164 local ETGs in Lauer et al. (2007a), despite the fact that the galaxy is more luminous than $\gtrsim 90\%$ of the sample ($M_{V}=-23.8\pm 0.1$, López-Cruz et al. 2014; see also Figure 12). It would also explain the large core size of Holm 15A.

Moreover, it could even provide a reason for Holm 15A’s large cusp radius (Figure 11e): In the merger simulation during the core-core re-merger $M_{BH}$ doubled while the core radius (described either by $r_{b}$ or $r_{\gamma}$) roughly tripled in size. This would suggest that in successive core scouring events the core grows faster then the central black hole. Similarly, for a sequence of five smaller core scourings due to minor mergers, the remnant also “outgrew” its black hole by a similar factor.

In the merger case, Holm 15A represents a dynamically very evolved galaxy that is possibly one merger generation ahead of cored galaxies like NGC4874 and NGC4889 at the center of the Coma cluster. As we showed in the previous subsection, Holm 15A’s high $M_{BH}/M_{Bu}$ ratio of $\sim 2\%$ might indicate that the galaxy’s progenitors had already formed at redshifts larger than 1 or 2 and/or that its progenitors were themselves BCGs. Abell 85 has one of the strongest cool-cores among X-ray bright clusters (Chen et al., 2007) and is strongly BCG dominated, with Bautz-Morgan morphological type I (Hudson et al., 2010) such that the central parts of the main cluster in fact might have been subject to a slightly accelerated evolution at some point in the past. Previous X-ray studies of Abell 85 had already suggested that the measured temperature and metallicity maps of the cluster were compatible with an intense merger history (e.g. Durret et al., 2005; McDonald et al., 2010).

5.5 Alternative formation scenario via AGN feedback?

Even though the merger scenario provides a consistent explanation for the central light profile shape of the galaxy, the orbital structure and how both are connected to the mass of the central black hole, we briefly discuss whether the interaction between an AGN and the surrounding stars could serve as an alternative core-formation scenario.

In recent simulations, AGN outflows have been observed to trigger fluctuations of the local gravitational potential which irreversibly transfer energy to the dark matter and stellar components (Teyssier et al., 2011; Martizzi et al., 2012, 2013; Choi et al., 2018). These simulations produced exponential light profiles, which resemble the cores of ETGs in the sense discussed in the introduction: the central surface brightness is low and the slope of the central surface brightness profile is shallow. In fact, based on the black-hole fundamental plane it has been argued that many black holes in the BCGs of cool-core clusters could be more massive than predicted by the classical black-hole scaling relations, and many would actually be expected to have masses $M_{\mathrm{BH}}>10^{10}\,M_{\odot}$ (Phipps et al., 2019; Hlavacek-Larrondo et al., 2012). We are still lacking numerical simulations that study in quantitative detail the effect of AGN feedback on the stellar light distribution and orbital structure. The information contained in the actual orbits of the stars might turn out to be crucial to distinguish between different core formation scenarios.

5.6 Dark matter halo and stellar mass-to-light ratio

Figure 15 shows the underlying stellar, dark matter and total enclosed mass and density profiles of our best-fit dynamical model of Holm 15A. Apart from the $20\%$ variation in $M_{BH}$, the quadrants of the galaxy produce a consistent overall mass and density profile.

Using simple stellar population models (Thomas et al., 2003; Maraston & Strömbäck, 2011) we find that Holm 15A has a marginally super-solar metallicity, $[Z/H]=0.08\pm 0.05$ and is strongly $\alpha$-enhanced $[\alpha/Fe]=0.25\pm 0.03$. Assuming a Kroupa stellar initial mass function (IMF) we find a stellar mass-to-light ratio of $\Upsilon_{SSP,Kroupa}=2.7\pm 0.30$ (i-band) using methods from either Maraston & Strömbäck (2011) or Conroy et al. (2017). The large $\sim 20\%$ uncertainty of this value is due to the difficulty of determining the age of the stars. Formally, our SSP models fitted stellar ages that exceed the age of the universe. The value of $\Upsilon_{SSP,Kroupa}$ and its uncertainty are derived from “manually” varying stellar ages between $10\text{\,}\mathrm{G}\mathrm{y}\mathrm{r}\mathrm{s}$ and $13.8\text{\,}\mathrm{G}\mathrm{y}\mathrm{r}\mathrm{s}$ while fixing elemental abundances.

Our dynamical mass-to-light ratio of $\Upsilon_{\star}=4.5\pm 0.19$ is roughly twice as large as the SSP ratio ($\Upsilon_{\star}/\Upsilon_{SSP,Kroupa}=1.7\pm 0.20$). This is a continuation of a growing trend among recent mass-to-light ratio measurements in massive ETGs from dynamics, lensing and spectroscopy often finding values larger than predicted by SSP models adopting a Kroupa stellar IMF, $\Upsilon_{\star}/\Upsilon_{SSP,Kroupa}\gtrsim 1.6$ (e.g. Treu et al., 2010; Auger et al., 2010; Thomas et al., 2011; Spiniello et al., 2011; Cappellari et al., 2012; Conroy & van Dokkum, 2012; Tortora et al., 2014; Conroy et al., 2017; Parikh et al., 2018; Alton et al., 2018). This offset is roughly consistent with a mass-to-light ratio implied by Salpeter-like IMF or might suggest that DM traces the stars. Our stellar-dynamical mass-to-light ratio is based on the assumption that all mass tracing the galaxy’s light profile belongs to the stars of the galaxy. In this case, when parameterizing the inner DM-halo as $\rho_{DM}\sim r^{-\eta}$, we find $\eta=0.45\pm 0.16$ out to roughly $50\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$. This is substantially shallower than predicted by numerical simulations of cold dark matter, $\eta\geq 1$ (e.g Navarro et al., 1996, 1997; Moore et al., 1998). Combined stellar kinematics and weak & strong lensing studies of local BCGs previously found $\rho_{DM}\sim r^{-0.5}$ on scales comparable to the effective radius (e.g Sand et al., 2004, 2008; Newman et al., 2013).

Within the core region the fraction of DM is $\lesssim 20\%$. However, under the assumption of a Kroupa IMF and that DM traces stars, the fraction of DM within the core region would be roughly $50\%$, while in the former scenario equality between the enclosed stellar and DM mass is reached only at $r_{eq}=33\pm$2.5\text{,}\mathrm{k}\mathrm{p}\mathrm{c} (The stellar mass density profile reaches equality with the DM density profile at $28\pm$0.10\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}). In both scenarios the mass density distribution of the stars in our best-fit model has a slope similar to that of the distribution of DM inside the core, $\rho_{total}\sim r^{-0.5}$.

We note that some massive galaxies seem consistent with a low-mass IMF (e.g. Thomas et al., 2016; Collier et al., 2018) and that some fine-tuning is required to consistently combine masses from multiple constraints like lensing, dynamics or spectroscopy (e.g. Newman et al., 2017). Dynamical and lensing constraints, in general, become model dependent when stars and DM trace each other closely (e.g. Thomas et al., 2011).

6 Summary and Conclusions

We have observed Holm 15A, the BCG of the cool-core galaxy cluster Abell 85, with MUSE. Our observations reveal a galaxy with little rotation ($v_{rot}<40\,\mathrm{km/s}$) and a nearly constant velocity dispersion of $\sigma=340\,\mathrm{km/s}$. Towards the center and towards large radii, the velocity dispersion increases slightly.

We use orbit-based, axisymmetric Schwarzschild models to analyse the dynamical structure of Holm 15A and compare them to recent high-resolution N-body simulations of mergers between ETG galaxies that host black holes. Our results indicate the following:

• •

Holm 15A hosts a $(4.0\pm 0.8)\times 10^{10}\$\mathrm{M}_{\odot}$$ SMBH at its center, the most massive black hole directly deteced via stellar dynamics so far. The black hole constitutes close to $2\%$ of the total stellar mass of the galaxy.

• •

Inside of the gravitational sphere of influence of the black hole, $r_{SOI}=3.8\pm$0.37\text{,}\mathrm{k}\mathrm{p}\mathrm{c}$$, the orbital distribution becomes increasingly tangentially anisotropic. However, the anisotropy inside the core is less tangential than in other big elliptical galaxies with depleted cores.

• •

The galaxy’s light profile and the observed mild orbital anisotropy both match remarkably well with predictions from N-body simulations of a merger between two elliptical galaxies that already had depleted cores.

• •

The SMBH is roughly 9 times larger than expected from the $M_{BH}-\sigma$ relation and 4 times larger than expected from the stellar mass of the galaxy, when compared to other cored ETGs. However, the offsets are smaller when compared to other BCGs.

• •

In core galaxies black hole masses scale inversely with the central stellar surface brightness $\mu_{0}$ and central stellar mass density $\Sigma_{0}$ - including in Holm 15A. We show this correlation here for the first time.

• •

Even in extreme instances of core formation like in Holm 15A, the core-specific relations $M_{BH}-\Sigma_{0}$, $M_{BH}-\mu_{0}$, $r_{SOI}-r_{\gamma}$, as well as the global galactic relation $M_{BH}-L_{K}$ still seem to hold. But the details of the light profile and orbital anisotropy contain valuable information about the specific formation path.

• •

Assuming that all the mass that follows the light is stellar, we infer a bottom-heavy IMF, $\Upsilon_{\ast}=4.5\pm 0.19$ ($i$-band), and the inner power-law slope of the DM-density distribution to be $\eta=0.45\pm 0.16$. Equality between enclosed stellar and DM mass is reached at $33\pm$2.5\text{,}\mathrm{k}\mathrm{p}\mathrm{c}$$. Assuming a Kroupa IMF, $\Upsilon_{SSP,Kroupa}=2.7\pm 0.3$, and DM tracing stars, we infer $\eta\sim 1$ outside of the core and a DM-fraction of nearly $50\%$ within the core.

We plan to extend our analysis of the galaxy to triaxial Schwarzschild models. This will allow us to investigate potential systematics related to symmetry assumptions in the modelling and related to possible substructure near the very center of the galaxy.

Our results suggest that the exact shape of the central light profile as well as the details of the distribution of stellar orbits in the center contain valuable information about the merging history of very massive galaxies. E.g., extreme instances of core formation could potentially lead to remnant surface-brightness profiles diverging from the typical core-Sérsic profiles of “classical” cored galaxies. Hydrodynamical cosmological simulations have also produced large stellar and dark-matter cores through AGN feedback. It will be interesting to compare the anisotropy profiles predicted by these simulations with measurements in observed galaxies. More extensive simulations are also required to investigate in detail the effect of core scouring under different initial conditions of the progenitor galaxies and on the DM halo.

The SMBH of Holm 15A is a candidate system for direct imaging of its sphere of influence. The photon ring radius is $\sqrt{27}GM_{BH}/c^{2}=2100\pm$410\text{,}\mathrm{A}\mathrm{U}. At redshift z = 0.055, this corresponds to an area spanning $18\pm$3.7\text{\,}\mu\mathrm{a}\mathrm{s} on the sky, only slightly smaller than the current minimum angular resolution of the Event Horizon Telescope, $25\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ (Event Horizon Telescope Collaboration et al., 2019).

Acknowledgement

We acknowledge the support by the DFG Cluster of Excellence "Origin and Structure of the Universe". The dynamical models have been done on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP) and we are grateful for the support by A. Krukau and F. Beaujean through the C2PAP. We are grateful to Hans Böhringer for valuable discussions and suggestions.

Appendix A Parametric analysis of Wendelstein Photometry

A.1 1D-Analysis of the Wendelstein image

The best-fit parameters of the various models we fit to the 1D $i$-band Wendelstein image of Holm 15A in Section 2.2 are shown in Table 2. In the table, the parameters of the different models are separated into components: Parameters of the Core-Sérsic function $I_{cS}$($r_{b}$, $n$, $\alpha$, $\gamma$, $n_{1}$, $r_{e,1}$) (see eq. 2 from Rusli et al. (2013a)), outer Sérsic function $I_{S}$($n_{2}$, $\mu_{e,2}$ and $r_{e,2}$) (cf. outer Sérsic components in Sérsic + Sérsic models from Kluge et al. (2019), eq. 11,12) and Nuker function $I_{N}$($r_{b}$, $n$, $\alpha$, $\beta$, $\gamma$) (see eq. 10 from Lauer et al. (2007b)).

A.2 2D-Analysis of the Wendelstein image

As described in Sec. 2.2, a detailed investigation of the 1D light profile of Holm 15A did not provide strong evidence for a break radius that separates the inner core from the rest of the galaxy. Here, we describe in detail our 2D fits to the $i$-band image using Imfit (Erwin, 2015). Our goal in performing these fits was to better understand the structure of the unusual core region of Holm 15A. In particular, whether or not a 2D analysis including the ellipticity structure of the galaxy would help in constraining the size of the galaxy’s core.

To have a fully independent analysis, we created ellipse fits to the Wendelstein image using the IRAF task ellipse (Carter, 1978; Jedrzejewski, 1987), complementary to our analysis in Section 2. The surface brightness profile and isophote shape measurements out to 250″ (see Figure 16) are fully compatible with the results from the other method (cf. Figures 1 and 2).

Beyond about 140″, the position angle twists by about 90°, and the ellipticity drops from $\sim 0.4$ to $\sim 0.2$ (see Figure 16). Meanwhile, the centers of the fitted ellipses begin varying by as much as $\sim 15\arcsec$. It is not clear how much of this represents a real change in the isophotes, e.g. if this is related to a transition to intra-cluster light, or how much is simply an artifact of the increasingly low S/N (signal-to-noise ratio). We therefore confine our 2D fitting to $a<140\arcsec$. In the intervall between $3\arcsec$ and at least $100\arcsec$, the position angle is remarkably stable, suggesting that Holm 15A might be close to rotational symmetry.

Towards the very center, the change in position angle implies that the isophotes start rotating but at the same time the galaxy becomes significantly rounder.

Because the region of the core is close to circular, the actual isophotes do not show any visible twists or distortions.

We find that fitting the image with an inner Sérsic function that is near-exponential in shape, with a Sérsic index of $n=0.99$, and an outer Sérsic component with $n=1.48$ results in a good fit to the Wendelstein data. The inner component is consistent with the Sérsic + Sérsic model listed in Kluge et al. (2019), though the outer Sérsic index is smaller. It is also smaller than for our core-Sérsic + Sérsic models from Section 2.2 (cf. Table 2). It is however consistent with models from Donzelli et al. (2011), who found that Holm 15A’s R-band light profile is well fit by the sum of two exponential functions (i.e. equiv. to a Sérsic + Sérsic model, with both $n\sim 1$). Similar results were obtained from the 2D analysis of a CFHT $r$-band image by Bonfini et al. (2015).

However, replacing the inner, exponential-like Sérsic component with a core-Sérsic component, did not significantly improve the quality of the fit relative to the core-less model. This reflects the radial trend of the observed light profile shown in the bottom panel of Figure 1 - The central light profile of Holm 15 is approximately exponential up to $\sim$25\text{,}\arcsec$$.

Nonetheless, there was still a distinct, bilobed excess in the residual image from the both the double Sérsic and the core-Sérsic + Sérsic fit, on a scale of $a\sim 4\arcsec$. We therefore experimented with adding additional components to the model. The best result was with the GaussianRing3D function of Imfit, which performs line-of-sight integration for an inclined ring with a Gaussian radial density and an exponential vertical density. The final result was a fit with central residuals which were almost completely lacking in any systematics (see Figure 17). The “ring” component has a semi-major axis of $4\farcs 1$, a position angle of 53° – almost perpendicular to the Sérsic components – and is intrinsically circular, viewed at an inclination of 68°. We also note that this may be consistent with the extra Gaussian-like Sérsic component (with $n=0.3$) – with a position angle of $\sim 55\arcdeg$ – which Bonfini et al. (2015) added to their 2D fits as a “corrective” component. We emphasize that this is a purely empirically chosen function which produces approximately the right excess light to minimize the residuals; it is not necessarily evidence for an actual inclined ring. The parameters of the best-fit 2D model are listed in Table 3.

In summary, while the 2D analysis provides somewhat more stable fit parameters, it confirms the results from Sec. 2.2, in particular the lack of a clear break radius. In the 2D analysis, we assume a spatially constant flattening for each individual component. This might imply that the components simply trace the structure of the ellipticity profile of the galaxy and this, in turn, could explain why the parameters of the 2D fits turn out more stable than in the 1D analysis. It is not clear at the moment how much or which physical information is encoded in the ellipticity profile of the galaxies. Likewise, it is not clear how physically significant the extra-light ring might be, which has a total luminosity comparable to the expected amount of stars ejected from the center by a SMBH binary, i.e. the extra light is of a similar order $\sim 0.5\times 10^{10}L_{\odot}$ as the missing light determined in Section 2.2. Comparing the distribution of stars in Holm 15A to those of other cored ETGs (cf. Figure 1) makes clear that Holm 15A is not only characterised by an extreme deficit of light in the inner core but also by an excess of light adjacent to the core. This light "excess", however, extends well beyond the extra-light ring (roughly out to $20\arcsec$).

Appendix B Stellar kinematics

B.1 Kinematics of Holm 15A compared to MASSIVE survey ETGs

To better understand Holm 15A’s place among other known massive ETGs we will compare it’s stellar kinematics to ETGs from the MASSIVE survey (Ma et al., 2014, and subsequent MASSIVE survey papers).

Characterizing Holm15 A’s velocity dispersion profile, $\sigma(r)$ (see Section 3.2) by fitting a combined power-law profile as suggested by Veale et al. (2018) in their study of the 90 ETGs of the MASSIVE survey, we find an inner logarithmic slope $\gamma_{inner}=-0.017\pm 0.007$ of the $\sigma$ profile at $\sim$2\text{,}\mathrm{k}\mathrm{p}\mathrm{c} and an outer logarithmic slope $\gamma_{outer}=0.029\pm 0.009$ at $\sim$20\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}. Roughly $90\%$ of BCGs in the MASSIVE survey have $\gamma_{inner}\leq 0$ and $\sim 60\%$ with $\gamma_{outer}\geq 0$. Moreover, for the eleven most massive BCGs in their sample with $M_{*}\sim 10^{12}$\mathrm{M}_{\odot}, $\gamma_{inner}\leq 0$ and $\gamma_{outer}\geq 0$ for all except one. The scatter in $\gamma_{inner}$ and $\gamma_{outer}$ between these eleven most massive BCGs is quite high, $\overline{\gamma}_{inner}=-0.040\pm 0.055$ and $\overline{\gamma}_{outer}=0.088\pm 0.084$. Nonetheless, statistically, their overall rather flat $\sigma(r)$ profiles are similar to the one in Holm 15A, even though the galaxy’s average velocity dispersion within one effective radius $\sigma_{e}\sim$340\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s} is slightly higher compared to these BCGs $\sim$300\text{,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$.

The parameter $h_{4}$, in our measured kinematic profile starts out at $\sim 0.07$ within $2\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$ and rises to $\gtrsim 0.1$ along the major axis towards the edges of the MUSE FOV. All 11 of the most massive MASSIVE BCGs share this trend of $h_{4}>0$ over their respective radial coverage and all but one have positive $h_{4}$ gradients towards larger radii. Similarly as with $\sigma$, average values for $h_{4}$ within $r_{e}$ are larger for Holm 15A, $h_{4,e}\sim 0.08$ than for those other BCGs where $h_{4,e}\lesssim 0.06$. Essentially all galaxies in the MASSIVE sample with $h_{4,e}>0.05$ (BCG or not) have within the central $2\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$ super-solar $[\alpha/Fe]>0.2$ and most galaxies with $h_{4,e}>0$ have $[Fe/H]\leq 0$ (Greene et al., 2019).

Using stellar population models of Lick indices (Thomas et al., 2003; Maraston & Strömbäck, 2011) we find abundance ratios in good agreement with these in Holm 15A: $[\alpha/Fe]=0.25\pm 0.03$ and $[Fe/H]=-0.011\pm 0.008$.

Overall, we find stellar kinematics in Holm 15A similar to those of other known massive ETGs. Indeed, from a stellar-kinematic point of view we find no indication that Holm 15A is anything other than a higher-mass extrapolation of known massive ETGs in the local universe, the vast majority of which is cored (e.g. Lauer et al., 2007a; Krajnović et al., 2013; Kormendy & Ho, 2013).

B.2 Non-parametric kinematics compared to Gauss-Hermite polynomials

We compare the non-parametric stellar kinematics we measured with our own code with those derived parametrically with pPXF. This is illustrated in Figure 18 for all bins of our FOV (i.e. LOSVDs from all quadrants). Both kinematic profiles are, for the purpose of illustration, parameterized via Gaussian times third to forth order Gauss-Hermite polynomials. As the distribution of differences in the right column of the figure show, both methods agree within their uncertainties.

Appendix C Non-parametric dynamical modeling: escape velocities

Here, we will briefly discuss the connection between the wings of the observed line-of-sight velocity distributions on the one side and the mass distribution and orbital structure on the other. Figure 19 shows an example of a non-parametric LOSVD measured near the center of Holm 15A together with the corresponding LOSVD from our best-fit dynamical model. We define the cutoff velocity $v_{0}$ of any LOSVD as the mean $v_{0}=(v_{0,+}+v_{0,-})/2$. If $v_{\mathrm{peak}}$ denotes the line-of-sight velocity at which the LOSVD peaks, then $v_{0,+}$ is the smallest zero of the LOSVD for $v_{\mathrm{los}}>v_{\mathrm{peak}}$ and $v_{0,-}$ is the absolute value of the largest zero of the LOSVD for $v_{\mathrm{los}}<v_{\mathrm{peak}}$, respectively. For Holm 15A this definition is sufficient since there is almost no detectable rotation and the LOSVDs are largely symmetric with respect to $v_{peak}$. For the LOSVD in Figure 19 we measure $v_{0}\sim$1375\text{,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$.

Figure 20 shows all the measured cutoff velocities $v_{0}$ from our MUSE observations together with the escape velocity curves $v_{\mathrm{esc}}(r)$ of the four best-fit models for the four quadrants of the galaxy. Here, we define $v_{\mathrm{esc}}$ relative to the maximum radius that is sampled by the orbit library. The uncertainties of the cutoff velocities are measured via the difference between values of $v_{0}$ determined from $LOSVD(v_{los})+\Delta LOSVD(v_{los})$ and $LOSVD(v_{los})-\Delta LOSVD(v_{los})$. Outside the core ($r\gtrsim 5\,\mathrm{kpc}$), the best-fit $v_{\mathrm{esc}}(r)$ curves follow closely the maximum observed cutoff velocities $v_{0}$. This is expected in a radially anisotropic system where a significant number of stars is populated on weakly bound, radially extended orbits. The less bound and the more radial the orbit is, the closer the orbital velocity gets to $v_{\mathrm{esc}}$. Indeed, outside the core region, our best-fit models become increasingly radially anisotropic (cf. Figure 14).

The situation changes towards the center of the galaxy, where the gravitational well is deepest. The observed cutoff velocities decrease at small radii, whereas the escape velocity necessarily increases. This can only be explained as an anisotropy effect: inside the sphere-of-influence of the central black hole (indicated by the vertical line), the orbit distribution becomes tangential (cf. figure 14). Since only stars on the most radial orbits can move with velocities up to the escape velocity and those stars are missing, the LOSVDs do not extend to $v_{\mathrm{esc}}$ anymore but vanish at smaller velocities.

The uncertainties in the observed cutoff velocities are large (due to the noise in the wings of the LOSVDs). This is indicated by the large scatter in values of $v_{0}$. However, the figure clearly demonstrates the importance of the information contained in the wings of the LOSVDs for both the gravitational potential as well as for the orbital structure.