Probing the Variation of the Inner-Surface-Brightness Profile of Nuclear Star Clusters on the Intermediate-Mass Black Hole Mass Measurements Using Mock Observations of ELT/MICADO and HARMONI

1. Introduction

There is a fascinating, but elusive, population of intermediate-mass black holes (IMBHs, $M_{\mathrm{BH}}\approx {10}^{2-5}$ ${\mathrm{M}}_{\odot }$) in low-mass galaxies [$M_{★}<{10}^{10}$ ${\mathrm{M}}_{\odot }$;][1,2,3,4,5,6], falling between the masses formed by single-massive stars today ($M_{\mathrm{star}}>100$ ${\mathrm{M}}_{\odot }$) and the supermassive black holes [SMBHs, $M_{\mathrm{BH}}\approx {10}^{6-10}$ ${\mathrm{M}}_{\odot }$;][7,8,9,10,11], which are often found at the centers of massive galaxies [$M_{★}\gtrsim {10}^{10-12}$ ${\mathrm{M}}_{\odot }$;][12,13]. The compelling evidence of IMBHs was the LIGO/Virgo discovery of an IMBH formed from a collision of two smaller black holes [BHs, $M_{\mathrm{BH}}={142}_{-16}^{+28}$ ${\mathrm{M}}_{\odot }$;][14]. However, current ground-based gravitational wave interferometers are limited to detecting BHs with $M_{\mathrm{BH}}<200$ ${\mathrm{M}}_{\odot }$. For $M_{\mathrm{BH}}\approx {10}^{3-5}$ ${\mathrm{M}}_{\odot }$, the search must be done with integral field spectroscopic (IFS) observations in nearby stellar systems [][hereafter N25][15] for stellar dynamics with the 39-meter European Extremely Large Telescope (ELT) and its High Angular Resolution Monolithic Optical and Near-infrared Integral (HARMONI) field spectrograph instrument [16].

It is also not yet clear if there’s a full spectrum of IMBH masses [17]. Nonetheless, understanding the abundance/small or rarity/massive of IMBHs contributes significantly to our cosmic inventory [18], providing fundamental insights into the structure of our universe. Previous studies have profound implications for understanding how SMBHs formed in terms of IMBH seeds in the early universe through several competitive channels like the dead remnants of the first-star generation [19], the gravitational runaway of densest stellar clusters [20], or the direct collapse of pristine gas within the massive halos [21,22]. They shed light on the connections with nuclear star clusters [NSCs;][23], the behavior of BH–galaxy scaling relations at the low-mass regime [24,25], the origins of ultracompact dwarf galaxies as former galaxy nuclei [26,27], and the detection rate predictions of gravitational wave of the Laser Interferometer Space Antenna [LISA;][28], which can detect longer-wavelength radiation than ground-based ones [29].

Another hypothesis for the massive BH seeds was the absence of direct observations of IMBHs with masses between those of stellar-mass BHs and SMBHs [30,31], leading to the speculation that IMBHs might not exist. However, this rarity is now considered a potential sample selection bias [23,24]. The sphere of gravitational influence (SOI, $R_{\mathrm{SOI}}=G{M}_{\mathrm{BH}}/{\sigma }_{★}^2$, where G is the gravitational constant and ${\sigma }_{★}$ is the central velocity dispersion of stars) around such IMBHs is typically too small to spatially resolve in most galaxies due to their further distance and small mass. For example, a IMBH with $M_{\mathrm{BH}}\approx {10}^5{M}_{\odot }$ and ${\sigma }_{★}\approx 30\mathrm{km}{\mathrm{s}}^{-1}$ within $3.5$ Mpc has $R_{\mathrm{SOI}}\approx 0_{.}^{″}03$, which is smaller than the diffraction limits of current state-of-the-art 8-10 meter class telescopes equipped with adaptive optics (AO) techniques, such as the Very Large Telescope (VLT) and Gemini (with a full width at half maximum (FWHM) of the point spread function (PSF), ${\mathrm{FWHM}}_{\mathrm{PSF}}\approx 0_{.}^{″}1$–$0_{.}^{″}2$), or Keck (≈$0_{.}^{″}05$).

However, there is now a growing number of IMBH candidates at the centers of dwarf galaxies, based on indirect estimates of their masses [32,33,34,35]. Chilingarian et al. [36] found 305 distant IMBH candidates with $M_{\mathrm{BH}}=(0.3-2)\times {10}^5$ M_⊙, ten of which were confirmed as active galactic nuclei (AGN) via coincident X-ray emission. Furthermore, a few IMBHs were directly detected in nearby dwarfs with bright NSCs [2,37,38,39,40,41,42,43,44,45] or in globular clusters [GCs;][46,47,48,49,50,51], although these detections are often at the limit of what 8-10 meter telescopes can achieve.

These observations are necessary the exploration of the relationships between $M_{\mathrm{BH}}$ and various galaxy properties, including the stellar mass [52], the luminosity [53], and the statistically velocity dispersion [54] of the central spherical bulge of stars. These properties, once measured in other galaxies, can serve as indicators of $M_{\mathrm{BH}}$. Therefore, firm detections and measurements of these currently elusive IMBHs in the local universe are crucial for robustly constraining the BH–galaxy scaling correlations across the mass scale. These strongly suggest a close connection between the growth of SMBHs and the formation and evolution of their host galaxies [55].

Given such importance of IMBHs and the current active races in searching for their existence like The Measuring Black Holes in Below Milky Way mass galaxies Project [The MBHBM_★;][44,56,57,58], N25 started a large survey for future searching for such IMBHs within NSCs within 10 Mpc by ultilising the unprecedented spatial and spectral resolutions and sensitivity of HARMONI. Later extended the sample to further distances up to the location of Virgo Cluster [$D\approx 16.5$ Mpc;][59] and Fornax Cluster [$D\approx 20$ Mpc;][60]. These works used the observed surface brightness (SB) from Hubble Space Telescope (HST) images, an assumed synthetic spectrum of their stellar population, Jeans Anisotropic Model [JAM1; 61,62] of the stellar dynamics, and the HARMONI simulator [HSIM2; 63] to create mock HARMONI observations. They analyzed the mock data cubes like real data, using JAM to infer the IMBH mass and its error in a Bayesian framework. While N25 simulations showed that the HARMONI instrument can clearly detect the existence of IMBH demographics in NSCs down to a mass range of $3\times {10}^3-3\times {10}^4$ ${\mathrm{M}}_{\odot }$ within 10 Mpc, such similar findings at the distances of 20 Mpc required the IMBH mass must increase up to ≈${10}^5$ ${\mathrm{M}}_{\odot }$ [60].

Nevertheless, when modeling the mock HARMONI observations, N25 utilized HST/Wide Field and Planetary Camera 2 (WFPC2) Planetary Camera (PC) imaging, which has a spatial resolution of $0_{.}^{″}05$, a factor of five coarser than their HARMONI IMBH survey spatial scale (10 milli-arcseconds, or mas). They used this data to model the galaxy stellar mass in JAM to predict the stellar movements under the influence of both the IMBH and the entire galaxy’s gravitational potentials. To achieve the required spatial scale for the stellar-mass model, they interpolated the one-dimensional (1D) SB to a scale of 11 times smaller, assuming the NSC component is well described by a core-Sérsic function with an inner power-law slope of $\gamma =0.1$. This assumption is weak because the current facilities cannot resolve the vicinity of IMBHs, so the actual light distribution of NSCs at this small scale is unknown but could be probed with the ELT. Thus, the assumption of $\gamma =0.1$ may not be universal and could vary widely, thus biasing the IMBH mass determinations.

In this work, we improve the stellar-mass model for galaxies using SimCADO3 to simulate mock observations from the Multi-AO Imaging Camera for Deep Observations (MICADO) imager [64,65]. We will then use these new two-dimensional (2D) high-spatial-resolution images, combined with the mock HARMONI kinematics from N25, to re-estimate the mass of IMBHs at the centers of nearby NSCs. This approach not only re-examines the feasibility of discovering IMBH populations (N25) for the upcoming large HARMONI IFS and MICADO imaging surveys, but also validates accurate IMBH measurements by exploring the effect of variations in the inner part of SB of NSCs.

We outline the use of SimCADO and describe how we employed it to model the I-band images expected to be observed with the MICADO imager for NGC 300 and NGC 3115 dw01, varying the inner-power-law slope ($\gamma$) of the core-Sérsic function used to represent the SB of these NSCs in Section 2 . Following this, we detail the mock kinematics measurements of these nuclei obtained from mock HARMONI IFS using HSIM from N25. We then explain the JAM, which we used to explore the effect of varying the inner SB of these NSCs on the IMBH mass measurements in Section 3 . Finally, we discuss and summarize our findings in Section 4 .

2. New Stellar Mass Model

2.1. MICADO Imager

MICADO is the first-generation ELT imager, designed to produce diffraction-limited imaging and long-slit spectroscopy at near-infrared (NIR) wavelengths. It operates with both a multi-conjugate laser guide star adaptive optics (MCAO) and a single-conjugate natural guide star adaptive optics (SCAO). MICADO can achieve a sensitivity comparable to that of the James Webb Space Telescope (JWST) and offers a resolution improved by a factor of six. The instrument can produce images covering a field-of-view (FoV) of $50.5\times 50.5{\mathrm{arcsec}}^2$ at 0.82–2.45 $\mu \mathrm{m}$ in wide-field mode, using an array of nine $4096\times 4096$-pixel detectors with a plate scale of $4\mathrm{mas}$. In zoom mode, with a FoV of $18\times 18{\mathrm{arcsec}}^2$, MICADO provides a plate scale of $1.5\mathrm{mas}$ in high-contrast imaging mode. The instrument can perform astrometry with $50\mathrm{mas}$ precision for brighter sources and enables coronagraphy using focal and pupil plane coronagraphs. Additionally, its observation modes include a long-slit spectrograph with a spectral resolution of $\lambda /\Delta \lambda \approx 2.0\times {10}^4$, covering two spectral ranges of 0.82–1.55 $\mu \mathrm{m}$ and 1.49–2.45 $\mu \mathrm{m}$.

MICADO can address a broad range of modern astrophysics. Its science drivers focus on themes such as galaxy evolution through observations of high-redshift galaxies and relic populations in local galaxies, the formation and evolution of galaxies in the early universe, and the dynamics of dense stellar systems. Additionally, it aims to explore the full mass spectrum of BHs in extragalactic galaxies and at the center of the Milky Way, including SMBHs in galaxy nuclei and IMBHs in dense stellar clusters (e.g., NSCs and GCs). Other key areas include studying the star formation history of galaxies through resolved stellar populations, characterizing exoplanets, observing planet formation and circumnuclear disks at small angular scales, and investigating various aspects of the Solar System.

2.2. SimCADO Package

SimCADO is a data simulation package [66] designed to generate realistic mock detector plane array read-out files for MICADO. This Python 3 routine models the incoming light’s optical path, representing elements in the optical train by taking into account the effects of the atmosphere, telescope, instrument, and detector [67]. It can simulate the MICADO imaging with 4 mas and 1.5 mas per pixel in the wavelength range of 0.7–2.5 $\mu \mathrm{m}$ and provide the users with raw data sets that will be similar to what MICADO will produce during a typical observing run.

SimCADO is also highly configurable in simulating various observational scenarios equipped with SCAO and MCAO and integrated with subsystem effects along the optical train (the performance of the derotator or atmospheric dispersion corrector). While the former effects are referred to spectral (1D of wavelength $\lambda$) and spatial (2D of ($x,y$) spatial position). The latter include the shifts, rotations, convolutions, and distortions, or the cross effect of both the spatial and spectral nature of the PSF for the incoming photons. Additionally, SimCADO includes some additionally realistic features for simulating mock imagery for the MICADO detector array like the world coordinate system, variable sky background, PSF variability over the FoV, instrumental distortion map, missing segments due to the mirror surface’s re-coat, and extra-terrestrial optical path elements because of Zodiacal light, galactic extinction, atmospheric extinction, and scattered moonlight.

In practice, SimCADO accepts inputs in FITS or ASCII formats, and outputs all data as standard FITS files, allowing flexibly with other programs already in use within the astronomical community.

2.3. Inner-Surface Brightness Profile Variation

2.3.1. Mock MICADO Image Simulations

We used SimCADO to simulate the I-band images for NGC 300 and NGC 3115 dw01 observed with MICADO in the wide-field mode. Here, we used the best-fit parameters of their NSCs’ SB constrained from the HST/WFPC2 F814W images as discussed in section 5 and summarised in table 4 of N25.

N25 adopt a combination of a core-Sérsic and a Sérsic model to fit the NSC and the extended component of these galaxies, respectively. The core-Sérsic model [68,69] parametrizes the projected SB of NSCs, rather than their intrinsic density, as follows:

$$I\left(r\right)=I^{\prime }{\left[1+{\left({\displaystyle \frac{R_b}{r}}\right)}^{\alpha }\right]}^{\gamma /\alpha }exp\left[-b_{n_1}{\left({\displaystyle \frac{r^{\alpha }+R_b^{\alpha }}{R_e^{\alpha }}}\right)}^{1/\left(\alpha n_1\right)}\right]$$

(1)

where

$$I^{\prime }=I_b{2}^{-\gamma /\alpha }exp\left[b_{n_1}{2}^{1/\left(\alpha n_1\right)}{\left({\displaystyle \frac{R_b}{R_e}}\right)}^{1/n_1}\right]$$

In this function, the Sérsic index ($n_1$) and the power-law slope ($\gamma$) control the shapes of the outer Sérsic part and the inner power-law regime, while the sharpness parameter ($\alpha$) determines the transition between these two regions. $R_b$ is the break radius at which the transition is happened, and thus $I_b$ is the intensity at $R_b$ (converted to SB ${\mu }_b$ in the legends of Figure 1 and Figure 2). Outside $R_b$, the profile follows a Sersic [70] profile with a projected half-light radius $R_e$, but it gradually transitions to a power-law SB $I\left(r\right)\propto r^{-\gamma }$ at smaller radii $r\ll R_b$. Finally, the parameter $b_{n_1}$ is given as $b_{n_1}\approx 2n_1-{\displaystyle \frac{1}{3}}+{\displaystyle \frac{4}{405n_1}}+{\displaystyle \frac{46}{25515n_1^2}}+\mathrm{O}\left({\displaystyle \frac{1}{n_1^3}}\right)$ [71].

On the other hand, the Sersic [70] model used to describe the extended disk of these galaxies has a form:

$$I\left(r\right)=I_{e}exp\left\{-b_{n_2}\left[{\left({\displaystyle \frac{r}{r_e}}\right)}^{1/n_2}-1\right]\right\}$$

(2)

where $n_2$ represents the Sérsic index, distinguishing it from $n_1$ of the core-Sérsic profile in eq. 1, and it relates to the coefficient $b_{n_2}\approx 2n_2-1/3$. Here, $r_e$ is the effective radius of the profile, and $I_e$ is the intensity at $r_e$, which is converted to SB ${\mu }_e$, as shown in the legends of Figure 1 and Figure 2.

The cusp in the inner slope ($\gamma$) of the core-Sérsic profile reflects our current lack of knowledge about the profiles of NSCs at small radii. These regions are inaccessible with current facilities but will become observable with the ELT. N25 adopted a power-law index of $\gamma =0.1$ in their extrapolated SB and simulations at these radii. In this work, we modelled the I-band MICADO images based on the best-fit parameters of these galaxies’ SB obtained from N25, but allowed $\gamma$ to vary: $\gamma =0,0.1,0.2,0.5,0.7$ for NGC 300 and $\gamma =0,0.1,0.2,0.4,0.7$ for NGC 3115 dw01.

We clarify the choice of $\gamma$ for both galaxies. The simulations with $\gamma =0$ serve to validate predictions where the inner SB approaches a constant value, and the velocity dispersion profile asymptotically becomes a constant positive value, rather than declining toward zero—even in the absence of an IMBH, as demonstrated in earlier numerical [72] and simulation work (N25). Simulations with $\gamma =0.1$ are to assess the consistency of the galaxy mass models generated from the mock MICADO I-band images in this study with those interpolated into the 4 mas regime, derived from the HST/WFPC2 F814W images. For $\gamma >0.1$, we explore how $\gamma$ variation affects the determination of IMBH mass through stellar dynamics when using MICADO data.

We enhanced the SimCADO package by adding a 2D core-Sérsic function (eq. 1), complementing the existing 2D Sérsic function (eq. 2). During our simulations, we first generated emission curves based on the internal spectral energy distribution (SED) for an elliptical galaxy (NGC 3115 dw01) and a spiral galaxy (NGC 300) in the I-band filter. We then adjusted the spectral dimension to longer wavelengths, corresponding to the redshifts of NGC 300 ($z=0.00048$) and NGC 3115 dw01 ($z=0.00227$), as reported on the NASA/IPAC Extragalactic Database (NED4). Although these redshifts are negligible, this step was included for accuracy.

For each galaxy, we created 2D models using the core-Sérsic profile for the NSC ($q=0.999$) and the Sérsic profile for the extended-disk component ($q=0.999$ for NGC 300 and elliptical with $q=0.9$ for NGC 3311 dw01), then combined these models into a single source. We also fixed certain observational parameters to mimic the realistic observations, including obs_dit = 900s (observing time per exposure) and n_dit = 4 (number of exposures).

2.3.2. PSF Determination from MICADO Point Source Images

Accurate characterization of the PSF is essential for deconvolving galaxy light distributions and performing reliable photometric and structural analyses. For MICADO, the PSF is expected to be complex and sensitive to both AO performance and observing conditions.

To model the PSF, we simulate MICADO observations of point sources (i.e., stars) using SimCADO, adopting the same instrumental setup and conditions as those used for the galaxy simulations. The resulting stellar images are analyzed with the Multi-Gaussian Expansion (MGE) method, implemented through the mge.fit_sectors routine from the MgeFit Python package5 [73]. The mge.fit_sectors algorithm models the 2D light distribution of each simulated star as a sum of co-axial Gaussian components. This simulated-MGE PSF is subsequently incorporated into the analysis of mock MICADO galaxy images to accurately account for instrumental and atmospheric broadening effects.

2.4. Brightness Profiles from MICADO Galaxy Images

We showed in Figure 1 for NGC 300 and Figure 2 for NGC 3115 dw01 the I-band MICADO images with varying the inner-power-law slope of the core-Sérsic, as well as their 1D SB extracted from the Image Reduction and Analysis Facility (IRAF) ellipse task [74]. The ellipse routine systematically integrated the flux of stars within concentric annuli, allowing for variations in position angles and ellipticities along the galaxy’s semi-major axis and deconvolving with the simulated PSF image (subsubsection 2.3.2). Subsequently, we converted the average flux within each annulus (measured in counts/s) into SB expressed in mag.arcsec⁻². The photometric zero point for the I band is zero_point_I = 29.492 mag obtained from the photometric zero point function within the SimCADO framework, which simulated the light from a point source (i.e., a stars) and calculated as zero_point_I = 2.5log₁₀(counts/t_exp) + input_mag; here, t_exp = 1s and input_mag = 10 mag.

To highlight the effects of varying $\gamma$ of the core-Sérsic profile on the SB at radii probed only by the ELT, we overlaid the HST/WFPC2 F184W IRAF elipse profiles and their best-fit core-Sérsic + Sérsic functions from N25 in the same figure for each galaxy. It is evident that at these unexplored radii, larger $\gamma$ significantly increase the light distribution of stars within the NSC towards its center. This change will indeed impact the determination of the IMBH mass (subsection 3.3). To validate our simulations, we compared the previous profile with $\gamma =0.1$ fitted by N25 and found these profiles consistent to each other in both galaxies. This consistency is also evident in the case of $\gamma =0$, where our simulations produced flat profiles toward the galaxy centers, as demonstrated in previous numerical work by Tremaine et al. [72] and simulations (fig. 4 and section 4.2 of N25).

We converted these MICADO I-band images to MGE models using the Python version of the MgeFit with the mge_fit_sectors_regularized routine. During the fits, we also convolved the images with the MICADO I-band PSF. We saved these MGE-light models of the two galaxies varying with the inner-power-law slope ($\gamma$) in Table 1 and illustrated in Figure 3. This Figure shows the agreement/disagreement between the data and the model as 2D contours at equivalent radii and contour levels for the mock MICADO simulations with five simulated $\gamma$ for each galaxy. These MGE models can be analytically de-projected with a specified axis ratio and inclination to reconstruct the three-dimensional (3D) distribution of the entire galaxy.

Finally, we converted the luminosity surface density into mass surface density by assuming a constant mass-to-light ratio ($M/L$). For NGC 300, we used $M/L_{\mathrm{F}814\mathrm{W},\mathrm{dyn}}\approx M/L_I\approx 0.6$ (${\mathrm{M}}_{\odot }$/${\mathrm{L}}_{\odot }$) as reported by Neumayer and Walcher [75], and for NGC 3115 dw01 we adopted $M/L_{\mathrm{F}814\mathrm{W},\mathrm{phot}}\approx M/L_I\approx 1.4$ (${\mathrm{M}}_{\odot }$/${\mathrm{L}}_{\odot }$) derived by Pechetti et al. [76].

3. Dynamical Model

3.1. Jeans Anisotropic Model (JAM)

We utilize the JAM model based on the solution of the Jeans equations, assuming axisymmetry with a cylindrical alignment of the velocity ellipsoi, which is axially symmetric around the vertical direction (${\sigma }_z\ne {\sigma }_r={\sigma }_{\varphi }$). This JAM-predicted LOSVD is calculated in combination with the new mass model simulated by SimCADO ( Section 2 ), is then compared with the mock kinematics generated by N25 through a Bayesian framework to infer the IMBH mass and its error.

This choice of cylindrical aligned JAM is to ensures consistency with the mock kinematics generated by N25, where they used the same model for their HARMONI IMBH sample of dwarf galaxies hosting bright NSCs with significant rotational characteristics. We implemented this by setting align=‘cyl’ in the jam_axi_proj procedure of the JamPy package [62].

3.2. Mock Stellar Kinematics

This work used the stellar kinematics measured from the mock data cubes produced from HSIM. The detailed descriptions of how these mock observations were produced in various spectral gratings of medium ($\lambda /\Delta \lambda \approx 7,100$; J, H, K) and high ($\lambda /\Delta \lambda \approx 17,400$; H-high, K-short, K-long) spectral resolution have presented in sections 6.2 and 6.3 of N25.

N25 showed these stellar kinematic measurements from their mock HARMONI datacubes simulated for three different $M_{\mathrm{BH}}$: = 0%$M_{\mathrm{NSC}}$, = 0.5%$M_{\mathrm{NSC}}$, and = 1%$M_{\mathrm{NSC}}$, in their figs 6 and 7 for the nuclei of NGC 300 and NGC 3115 dw01, respectively. These LOSVD maps include the rotation (V), velocity dispersion (${\sigma }_{★}$), and root-mean-squared velocity ($V_{\mathrm{rms}}=\sqrt{V^2+{\sigma }_{★}^2}$). We note that NGC 300 and NGC 3115 dw01 host at their centres an NSC with a mass of $M_{\mathrm{NSC}}$ $\approx {10}^6$ ${\mathrm{M}}_{\odot }$ and $M_{\mathrm{NSC}}$ $\approx 7\times {10}^6$ ${\mathrm{M}}_{\odot }$, respectively.

These stellar kinematic maps (i.e., with an assumed $\gamma =0.1$) show distinct properties of the presence or absence of an IMBH, exhibiting a distinctive central drops in ${\sigma }_{★}$ and $V_{\mathrm{rms}}$ (if nuclei are less rotating) maps when either no BH ($M_{\mathrm{BH}}=0$ ${\mathrm{M}}_{\odot }$) is present or when the $M_{\mathrm{BH}}$ are insufficient (e.g., $M_{\mathrm{BH}}=5\times {10}^3$ ${\mathrm{M}}_{\odot }$ for NGC 300 and $M_{\mathrm{BH}}=3.5\times {10}^4$ ${\mathrm{M}}_{\odot }$ for NGC 3115 dw01; or $M_{\mathrm{BH}}$ = 0.5%$M_{\mathrm{NSC}}$). This central-${\sigma }_{★}$ drop is a common feature in the predicted stellar kinematics of galaxies when no or very small IMBHs are assumed, and it holds true for a range of assumed anisotropic profiles and density profiles with $\gamma >0$ regardless of the Sérsic index (Tremaine et al. [72]; N25). This well-established feature of realistic galaxy models initially appeared to contradict early observations [77]. However, the apparent discrepancy was resolved when it was discovered that all massive galaxies contain SMBHs. Conversely, models featuring relatively large IMBHs generate centrally rising peaks in the ${\sigma }_{★}$ map at several central spaxels in the Keplerian maner ${\sigma }_{★}\propto r^{-1/2}$ (for $0<\gamma <2$) irrespective of the profile slope and Sérsic index, which correspond to the BH’s SOI, showing a clearly upturned trending as $M_{\mathrm{BH}}$ increases. This observation aligns with the general expectation that central ${\sigma }_{★}$ should increase in a Keplerian manner within a tiny region near the central IMBH, where its gravitational potential dominates over all others.

Given the highly consistent among N25’s mock stellar kinematics measurements across six spectral gratings and two spectral resolutions and the purpose of reducing the calculations, we adopt the stellar kinematics produced from the H-high grating to test the variation effect of the inner SB of NSCs simulated with MICADO on IMBH mass measurements only.

3.3. Improve IMBH Mass Measurements with MICADO

We performed fitting of the JAM to mock H-high kinematics ($V_{\mathrm{rms}}$), aimed at testing the variation effect of the inner-power law’s slope ($\gamma$) of the core-Sérsic function of NSCs simulated with MICADO on IMBH mass determinations and at constraining the model’s dynamical parameters. Here, we replaced the interpolated HST mass model toward the central 4-mas region produced by N25 by the 2D MICADO images simulated from the SimCADO with varying a range value of $\gamma$ (subsection 2.3).

The JAM is characterized with a central IMBH mass ($M_{\mathrm{BH}}$), the stellar orbitals represented by the anisotropy (${\beta }_z$), the galaxy morphological inclination angle (i), and the constant stellar-mass-scaling factor ($M/L$_I). We adopted a logarithmic scale for $M_{\mathrm{BH}}$ to ensure a more Gaussian sampling across multiple orders of magnitude, while other three parameters are sampled in linear scales. During the fitting process, JAM generated kinematic models that could be compared to their corresponding simulated values ($V_{\mathrm{rms}}$) while accounting for the HARMONI LTAO PSF with a ${\mathrm{FWHM}}_{\mathrm{PSF}}\approx 12$ mas.

To obtain the best-fit JAM, we utilized a Markov chain Monte Carlo (MCMC) simulation to explore the parameter space of $M_{\mathrm{BH}}$, $M/L$_I, ${\beta }_z$, and i, fitting the mock kinematic data and determining the best-fit parameters and uncertainties (statistic and kinematic measurement), using the adaptive Metropolis algorithm [78] within the Bayesian framework [adamet6; 79]. Our MCMC chains comprised of $3\times {10}^4$ iterations, with the initial 20% of iterations excluded as a burn-in phase. The full probability distribution function (PDF) was constructed from the remaining 80% of calculations. The best-fit parameters are the highest likelihood of the PDF, while the uncertainties for all four parameters were determined at 1$\sigma$ and 3$\sigma$ confidence levels, representing the ranges (16–84)% and (0.14–99.86)% of the PDF, respectively. We fixed the search ranges for four model parameters as follows:

• $lg(M_{\mathrm{BH}}/$ ${\mathrm{M}}_{\odot }$): from 0 to 6 (or $M_{\mathrm{BH}}$: from 0 to ${10}^6$ ${\mathrm{M}}_{\odot }$)

• $M/L$_I: from 0.1 to 3 (${\mathrm{M}}_{\odot }$/${\mathrm{L}}_{\odot }$)

• ${\beta }_z$: from −1.0 to 0.99

• i: from 5° to 90°

Initial guesses for the JAM parameters were fixed with ($M/L$_I, ${\beta }_z$, i) = (0.6, 0, 42) for NGC 300, while those = (1.4, 0, 42) for NGC 3115 dw01. The initial guesses for $M_{\mathrm{BH}}$ as those used in the simulated data cubes with $M_{\mathrm{BH}}=0,5\times {10}^3,{10}^4$ ${\mathrm{M}}_{\odot }$ for NGC 300 and $M_{\mathrm{BH}}=0,3.5\times {10}^4,7\times {10}^4$ ${\mathrm{M}}_{\odot }$ for NGC 3115 dw01.

In Figure 4 for NGC 300 and Figure 5 for NGC 3115d w01, the best-fit JAM parameters and their associated uncertainties are shown. These parameters describe the derived stellar kinematics for three different $M_{\mathrm{BH}}$ and five inner power-law slopes $\gamma$. The figures feature 2D scatter plots for each parameter, with colored points indicating their likelihood. White represents the maximum likelihood, while black signifies a confidence level smaller than 3$\sigma$. Additionally, histograms are provided to display the 1D distributions for each parameter. These 1D distributions were used to determine the best-fit values and their corresponding uncertainties, which are listed in Table 2 for NGC 300 and Table 3 for NGC 3115 dw01. The insert plots at the top-right corner of each PDF directly compare the best-fit JAM model and the mock $V_{\mathrm{rms}}$ on the same velocity scale.

One notable feature of the stellar velocity dispersion produced by the best-fit JAM is clearly illustrated in Figure 4 for NGC 300, in the case with no BH. However, this pattern is less evident in the corresponding models for NGC 3115 dw01, as shown in Figure 5. As the $\gamma$ value increases, the central velocity dispersion drop becomes more pronounced, with a larger extent and greater depth, requiring a larger $M_{\mathrm{BH}}$ for JAM to reproduce the mock $V_{\mathrm{rms}}$. An increase in $\gamma$ also raises the central stellar mass density, and JAM utilizes the $M/L$ parameter to account for the degeneracy between stellar mass density and $M_{\mathrm{BH}}$. We explore these complexities to explain our results in detail below.

For the tests with $\gamma =0$ (flat core; see the top-row panels of Figure 4 for NGC 300 and Figure 5 for NGC 3115 dw01), JAM cannot replicate the central drop in velocity dispersion (and subsequently in the $V\mathrm{rms}$ maps) for cases with input $M_{\mathrm{BH}}=0$ ${\mathrm{M}}_{\odot }$. This is because the projected dispersion predicted by JAM remains flat towards the centers when $\gamma =0$, e.g., see fig. 4 of N25 and Tremaine et al. [72], leading to biased lower $M_{\mathrm{BH}}$ estimates and resulting in an upper limit (NGC 300) or lower $M_{\mathrm{BH}}$ estimate (NGC 3115 dw01) for the mock kinematic cases with input $M_{\mathrm{BH}}=0.5\%$ or $1\%$ $M_{\mathrm{NSC}}$. In these scenarios, there is no cusp in the SB, which is dominated by stars with radii around the break radius, and the density follows $\rho \propto r^{-1}$. Consequently, when there is no BH, the depth of the central potential well is finite, and the projected dispersion remains asymptotically constant as $r\to 0$, dominated by stars inside the break radius. If a BH is present, the projected dispersion increases at a rate between ${\sigma }_{★}\propto r^{-1}$ and $r^{-0.5}$.

In all tests with the weak cusp models ($0<\gamma \le 0.2$), we successfully reproduced the results from N25 with $\gamma =0.1$, achieving uncertainties within 30% for $M/L$_I and 10% for $M_{\mathrm{BH}}$. This outcome was expected, as the MICADO images were modeled using the HST SB constrained in N25. Additionally, we validated the numerical prediction for $\gamma =0.2$ in NGC 300, as reported by N25, by replicating it with our mock MICADO simulations. This further confirms the accuracy of same $\gamma$ simulations for NGC 3115 dw01 in this work.

However, in the case of NGC 300, there is an opposite trend in $M/L$_I and $M_{\mathrm{BH}}$ compared to their trends in N25. Specifically, our JAM with the MICADO mass models provided lower $M_{\mathrm{BH}}$ but higher $M/L$_I than those in N25 for the same mock kinematics. This discrepancy is due to the MICADO images providing more light (and thus stellar mass) in the radius range of 0$0_{.}^{″}2$–$3_{.}^{″}$0 than the best-fit surface brightness profiles constrained from HST images, as seen in Figure 1. In the case of NGC 3115 dw01, the trend of $M/L$_I and $M_{\mathrm{BH}}$ is entirely consistent with the findings from N25, due to the consistency between its MICADO and HST images (see Figure 2). Here, the dynamical model compensates for the central velocity dispersion drop primarily through the negative covariance between $M_{\mathrm{BH}}$ and $M/L$_I. The increase in central stellar mass density remains either lower than or comparable to the $M_{\mathrm{BH}}$, thus having a negligible impact on $M_{\mathrm{BH}}$ determinations. Specifically, for NGC 300, the difference in central stellar mass is $M_{★,\gamma =0.2}-M_{★,\gamma =0}\approx 1.4\times {10}^4$ ${\mathrm{M}}_{\odot }$, and for NGC 3115 dw01, it is $M_{★,\gamma =0.2}-M_{★,\gamma =0}\approx 7.4\times {10}^4$ ${\mathrm{M}}_{\odot }$.

The 2D PDFs show the 3$\sigma$ confidence level “banana shape” of the negative covariance between $M_{\mathrm{BH}}$ and $M/L$_I. This covariance occurs due to the interplay between the gravitational potentials of central BHs and their host galaxies probed by both MICADO and HARMONI, where larger $M_{\mathrm{BH}}$ correspond to smaller $M/L$_I, and vice versa.

For the immediate cusp models ($0.2<\gamma \le 0.5$), we considered $\gamma =0.5$ for NGC 300 and $\gamma =0.4$ for NGC 3115 dw01. In these cases, the stellar mass density difference between these $\gamma$ profiles and the profile with $\gamma =0$ becomes comparable to, or even exceeds, the BH mass within $0_{.}^{″}1$ of the core-Sérsic profile of the NSCs. Specifically, $M_{★,\gamma =0.5}-M_{★,\gamma =0}\approx 4.25\times {10}^4$ ${\mathrm{M}}_{\odot }$ for NGC 300 (Figure 1) and $M_{★,\gamma =0.4}-M_{★,\gamma =0}\approx 7.73\times {10}^5$ ${\mathrm{M}}_{\odot }$ for NGC 3115 dw01 (Figure 2), leading JAM predicts lower $M_{\mathrm{BH}}$ or fail to detect it (see Figure 4 and Figure 5).

For the strong cusp models ($\gamma >0.5$). Here, we considered $\gamma =0.7$ for both NGC 300 and NGC 3115 dw01. In these models, the central stellar mass density of the NSC significantly surpasses the $M_{\mathrm{BH}}$. For NGC 300, the increase in stellar mass is $M_{★,\gamma =0.7}-M_{★,\gamma =0}\approx 2.0\times {10}^6$ ${\mathrm{M}}_{\odot }$ (Figure 1), and for NGC 3115 dw01, it is $M_{★,\gamma =0.7}-M_{★,\gamma =0}\approx 1.55\times {10}^6$ ${\mathrm{M}}_{\odot }$ (Figure 2). As a result, the $M_{\mathrm{BH}}$ is completely diminished.

In addition, it should also be noted that the variation of $\gamma$ does not impact the constraints on the inclination (i), which remains weakly constrained in our dynamical models due to the round shape of NSCs. Similarly, the orbital anisotropy (${\beta }_z$) is only slightly negative for NGC 300 and positive for NGC 3115 dw01 but close to zero, the assumed value in the N25 HSIM simulations, indicating minimal effect.

4. Conclusions

We have presented ELT/MICADO simulations for novel galaxy stellar model for measuring IMBH mass measurements using N25’s mock HARMONI kinematics. We summarize our conclusions as follows:

(i) We enhanced the SimCADO software, used for simulating MICADO images on the ELT, by incorporating a core-Sérsic function to model the SB of NSCs.

(ii) Our mock MICADO I-band images for NGC 300 and NGC 3115 dw01, assuming an inner power-law slope of $\gamma =0.1$ for the core-Sérsic profile, are fully consistent with the constraints from N25. This holds true at both the larger scales observed with HST and the interpolated 4 mas scale. Given this consistency for $\gamma =0.1$, we used it as a reference for further simulations with slopes of $\gamma =0,0.2,0.5,0.7$ for NGC 300, and $\gamma =0,0.2,0.4,0.7$ for NGC 3115 dw01.

(iii) These mock MICADO images (4-mas pixel size) for both galaxies can be achieved with a one-hour on-source exposure time. Since their SB remain consistent with the existing HST observational scale at large radii, increasing the exposure time to four hours did not result in significant image quality improvement. Therefore, we set one hour as the maximum on-source exposure time for MICADO observations in the HARMONI IMBH sample.

(iv) We used mock MICADO images to reconstruct the stellar mass models for both galaxies. These models, combined with mock HARMONI kinematics from N25, served as inputs for JAM to re-estimate the masses of their central IMBHs and measurement uncertainties. The best-fit recovered JAM, based on the assumed $\gamma =0$ mass profile (e.g., flat core models) regardless of Sérsic index, showed a constant positive value in the velocity dispersion toward the galaxy centers, rather than the expected central drop for cases with zero or small IMBHs, as predicted by Tremaine et al. [72] and N25, leading to an underestimation of the IMBH mass.

(v) Our same tests with assumed $\gamma >0$ mass profiles, and increasing $\gamma$ for the NSCs, exaggerate the central drop in both the size and depth of the velocity dispersion at a few central spaxels, where the stellar kinematics are dominated by the gravitational potential of the IMBH. This feature was also shown in fig. 4 and discussed in section 4.2 of N25, though based on numerical predictions. Here, we confirmed the $\gamma$-dependence of stellar kinematics by simulating realistic photometric observations using MICADO, in combination with mock HSIM IFS HARMONI kinematics and JAM. While increasing $\gamma$ causes significant drops in central velocity dispersion, it also raises the stellar mass density in that region (within $0_{.}^{″}1$) compared to that of the flat core models, which is accounted for and balanced in the dynamical models.

• Weak cusp ($0<\gamma \le 0.2$): while the central stellar mass density increases, it remains lower than or comparable to the $M_{\mathrm{BH}}$. Thus, JAM compensates for the drop in central velocity dispersion primarily through the negative covariance between $M_{\mathrm{BH}}$ and $M/L$.

• Immediate cusp ($0.2<\gamma \le 0.5$): the central stellar mass density increases significantly and even surpasses the $M_{\mathrm{BH}}$. JAM thus tend to underestimate the IMBH mass or none detection.

• Strong cusp ($\gamma >0.5$): the central stellar mass density becomes dominant over the $M_{\mathrm{BH}}$. As a result, the BH $M_{\mathrm{BH}}$ is entirely suppressed in JAM.

Software: 

Python 3.12: [80], Matplotlib 3.6: [81], NumPy 1.22: [82], SciPy 1.3: [83], photutils 0.7: [84], AstroPy 5.1 [85], AdaMet 2.0 [79], JamPy 7.2 [62], MgeFit 5.0 [73], and SimCADO [67].