© The Authors 2025

1. Introduction

As a compact, powerful engine at a galaxy’s center, supermassive black holes (SMBHs) play a crucial role in driving the evolution of their host galaxies (Kormendy & Ho 2013). Accurately measuring SMBH masses (M_BH) is essential for constraining galaxy–black hole correlations and distinguishing their various coevolution pathways (Graham & Sahu 2023).

The M_BH estimate for NGC 4736 (M 94) has long been debated and varies significantly depending on the method and assumed distance (4.3–6 Mpc). The Graham (2023a)M_BH–σ relation predicts M_BH ≈ (1.05 ± 0.64)×10⁷ M_⊙ based on a bulge stellar velocity dispersion of σ = 110 ± 5 km s⁻¹ (Barth et al. 2002). In contrast, Hubble Space Telescope (HST) observations of the broad Hα emission component yielded M_BH ≈ 3 × 10⁴ M_⊙ (Greene & Ho 2005). Mid-infrared correlations based on [Ne V] and [O IV] emission lines suggested M_BH =(2.3 ± 0.4)×10⁵ M_⊙ and (1.7 ± 0.5)×10⁵ M_⊙ (Dasyra et al. 2008), respectively. Furthermore, the X-ray luminosity (2–10 keV band) of the point-like source spatially coincident with the nucleus is log₁₀L_X = 38.48 erg s⁻¹ (Eracleous et al. 2002; Pellegrini et al. 2002; van Oers et al. 2017). Together with the 5 GHz radio luminosity of the unresolved radio nucleus (log₁₀L_R = 35.51 erg s⁻¹; van Oers et al. 2017), the Merloni et al. (2003) fundamental plane (FP) relation estimates M_BH$=8.2_{-1.8}^{+2.1}\times {10}^5$ M_⊙. However, alternative versions of the FP relation more suitable to the low-luminosity active galactic nucleus (LLAGN) regime predict higher masses: M_BH$=1.1_{-0.7}^{+1.2}\times {10}^7$ M_⊙ according to Plotkin et al. (2012), and M_BH$=1.6_{-1.0}^{+2.4}\times {10}^7$ M_⊙ according to Gültekin et al. (2019).

This work presents a dynamical stellar-based M_BH measurement for NGC 4736 using data observed by the James Webb Space Telescope (JWST) Near-Infrared Spectrograph (NIRSpec) G235H/F170LP grating. These high-quality integral field unit (IFU) data enable precise 2D stellar kinematics within a 3″ × 3″ field of view (FoV), resolving the longstanding debate on its M_BH. Our result also demonstrates that JWST can probe mass ranges previously thought to be beyond its reach.

2. Observations and results

2.1. JWST IFU

NGC 4736 was observed with NIRSpec G235H F170LP (PI: Anil C. Seth; PID: 02016) on February 2, 2023. The data consist of eight exposures, each 218.8 seconds long, for a total of 1750.4 seconds on source. The observation is centered on the galaxy’s nucleus, utilized a high-spectral-resolution configuration (R ∼ 2700), and covered the wavelength range 1.66–3.17 μm. The NRSIRS2RAPID readout pattern was used with 14 groups per integration, one integration per exposure, and a four-point medium-cycling dither. The data were reduced following the JWST calibration pipeline for NIRSpec IFU using the STScI pipeline v1.8.2 with the Calibration Reference Data System context 1063¹ (Perna et al. 2023). The final cube was produced with a pixel scale of $0\stackrel{″}{.}1$, which is 1.5 times smaller than the point spread function (PSF) full width at half maximum (FWHM) of $0\stackrel{″}{.}15$ at 2.3 μm (Fig. 7 of D’Eugenio et al. 2024). The intensity map of the reduced cube is shown in panel A of Fig. 1.

Fig. 1. Stellar kinematic measurements extracted from the JWST/NIRSpec G235H/F170LP IFU. Panel A: Logarithmic integrated intensity collapsed along the spectral dimension (excluding the detector gap at ∼2.41–2.49 μm). Panel B: PPXF fit for the central VORONOI bin. Part of the central bin spectrum, which contains the stellar absorption features at ∼2.3 μm (CO band heads), is plotted in black, and the best-fit empirical X-shooter stellar library (XSL; Appendix A.1) is overlaid in red. The residuals (data-model) are illustrated in green. Two vertical gray lines indicate the spectral region where we actually fit the model templates to the spectra of all bins. Panels C, D, and E: Rotation (V), velocity dispersion (σ), and root-mean-squared velocity ($V_{\mathrm{rms}}=\sqrt{V^2+{\sigma }^2}$). White contours show the intensity, which decreases by 0.5 magnitudes from the center outward. We also performed the same measurements using the PHOENIX high-spectral-resolution synthetic spectral template (Appendix A.2), which are shown in Fig. A.1.

We measured the stellar kinematics from the final cube, using part of the CO band heads’ stellar absorptions (2.20–2.36 μm). We first applied the adaptive VORONOI² binning method (Cappellari & Copin 2003) to spatially combine fluxes of many spaxels into a single bin, reaching the targeted-S/N = 50 per spectral pixel. Second, we performed log-re-binning along the spectral dimension for each binned spectrum with a velscale = 50 km s⁻¹ per pixel calculated using Eq. (8) of Cappellari (2017) and then applied the PENALIZED P IXEL-FITTING (PPXF) algorithm³ (Cappellari 2023) to fit the stellar spectral libraries (Appendix A) to these binned spectra; however, see Appendix B for our derivation and use of the line-spread function (LSF).

In low-ionization nuclear emission-line region (LINER) galaxies such as NGC 4736, warm dust emission is minimal in the near-infrared (NIR; see Fig. 7 of Ho 2008). Further evidence of this is provided by recent JWST/Mid-Infrared Instrument (MIRI) spectra, which show that the dust emission tends to be negligible toward shorter wavelengths (i.e., < 5 μm; Goold et al. 2024) in the nuclear regions of two nearby LLAGNs, Sombrero and NGC 1052. We fit the kinematic parameters V (rotation), σ (dispersion), h₃ (skewness), and h₄ (kurtosis) by setting moments = 4, and excluded the additive and multiplicative polynomials by setting degree = −1 and mdegree = 0 (Cappellari 2023). We accounted for the instrumental broadening of the NIRSpec spectra by convolving the stellar templates (Appendix A) with the instrumentally different dispersion between the data and the templates. The best-fit template superimposed on sections of the central bin’s spectrum is shown in panel B of Fig. 1.

The rotation (V), velocity dispersion (σ), and root-mean-squared velocity ($V_{\mathrm{rms}}=\sqrt{V^2+{\sigma }^2}$) are displayed in panels C, D, and E of Fig. 1, respectively. The JWST nuclear σ measurement is in full agreement with that derived from HST Space Telescope Imaging Spectrograph (STIS) data at a scale of $0\stackrel{″}{.}1$ (Pagotto et al. 2019), indicating the presence of a central velocity dispersion peak. This is the second clear stellar kinematic signature of a central SMBH, after the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI) observations of NGC 3640 (Thater et al. 2019). It is characterized by a strongly resolved central drop toward the galaxy center and a sharp rise within the innermost spaxels in both the σ and V_rms maps at the scale of the SMBH’s sphere of influence (SOI). The nucleus of NGC 4736 has low rotational characteristics, V/σ ≲ 0.1, across the NIRSpec FoV, but this is likely mainly due to the nearly face-on view of the galaxy disk. The detection in HST images of a strong, point-like UV source ≈2$\stackrel{″}{.}$5 northwest of the nucleus led to the intriguing suggestion that the source could be a secondary nuclear black hole, leftover from a galactic merger (Maoz et al. 2005). However, our stellar kinematic results do not reveal any signatures at that location. We validated our kinematic measurements by comparing them with those obtained by the PPXF fit to 1.67–1.83 μm (using [Mg I] λ1.711–1.737 μm) and 2.13–2.30 μm (using S Iλ2.188, 2.183, 2.179 μm). Our experiments demonstrate that the differences in kinematic measurements among these three fitting approaches are ≲3%.

2.2. HST imaging

We used the WFPC2 F814W (pixel scale of 0$\stackrel{″}{.}$1) observations of NGC 4736 taken on July 7, 2002 (Project ID 9042; PI: Smartt) to constrain the galaxy’s stellar mass model. These data include two exposures totaling 460 seconds. We also tested the systematics arising from the use of multiband HST imaging (see Appendix D).

The stellar-light distribution was approximated using a multi-Gaussian expansion (MGE) model⁴ ((Emsellem et al. 1994 and the PYTHON version of the mge_fit_sectors_regularized procedure (Cappellari 2002) and assuming a distance of D = 4.66 Mpc to NGC 4736 (i.e., the distance to the tip of the red giant branch; Karachentsev et al. 2003). We masked the HST image to exclude regions obscured by dust as described in Nguyen et al. (2020). Thus, we decomposed the PSF of the HST/WFPC2 F814W filter, which was computed using the TINYTIM package⁵ (Krist et al. 2011), into an MGE. This PSF MGE was then used as an input to the MGE fit of the F814W masked image. This masked HST image was fit with a sum of 2D Gaussians, which were convolved with the PSF MGE, and were then analytically deprojected with a given free inclination (i) to obtain a 3D axisymmetric light distribution. We list the parameters of each spatially deconvolved Gaussian of the MGE model in Table 1 and compare them with the F814W image in the left panel of Fig. C.1.

3. SMBH mass with the Jeans anisotropic model

Using a constant M/L_F814W for our stellar-mass model, the initial optimization yielded the best-fit Jeans anisotropic model (JAM; Appendix E) with M_BH = (1.7 ± 0.18)×10⁷ M_⊙ and M/L_F814W = 8.30 ± 0.05 M_⊙/L_⊙. While this model fits the data well within a radius of $0\stackrel{″}{.}75$, discrepancies emerge beyond this radius (left panel of Fig. C.2). We addressed this issue by adjusting M/L_F814W = 8.4 M_⊙/L_⊙, improving the agreement between the model and the data from $0\stackrel{″}{.}75$ to $2\stackrel{″}{.}0$. This correction reduced the M_BH to 1.55 × 10⁷ M_⊙ and resolved the discrepancies (see Fig. 2).

Fig. 2. Comparisons of the stellar kinematics measured from the NIRSpec G235H/F170LP observations and a range of JAMs. The black line shows the model with a constant M/LF814W, and colored lines the models that incorporate a radially varying M/LF814W plus 1σ error. The dashed vertical line marks the extent of the SMBH’s SOI. The data profile was obtained by averaging over circular annuli, and the associated errors are defined as the median of the Vrms uncertainties calculated by PPXF (i.e., ΔV and Δσ) of all bins within the circular annuli.

In the second optimization of the stellar mass model, we considered a M/L_F814W varying with radius. Specifically, we assumed a M/L profile that is constant at radii < 0$\stackrel{″}{.}$75 (${M/L}_{0\stackrel{″}{.}75}$) and > 2$\stackrel{″}{.}$0 (${M/L}_{2\stackrel{″}{.}0}$) and varies linearly between the two radii. The value of (${M/L}_{0\stackrel{″}{.}75}$) was fixed to 8.3 M_⊙/L_⊙ based on the best-fitting JAM model, while (${M/L}_{2\stackrel{″}{.}0}$) was left as a free fit parameter. Thus, the M/L(σ) for each Gaussian with dispersion sigma is

$$M/L(\sigma )={(M/L)}_{0\stackrel{″}{.}75}+(\sigma -0\stackrel{″}{.}75)\times \left[\frac{{(M/L)}_{2\stackrel{″}{.}0}-{(M/L)}_{0\stackrel{″}{.}75}}{1\stackrel{″}{.}25}\right].$$

We obtained the best-fit constraint with M_BH = (1.60 ± 0.16)×10⁷ M_⊙, ${M/L}_{2\stackrel{″}{.}0}=8.39\pm 0.04$ M_⊙/L_⊙, orbital anisotropy β_r = −0.52 ± 0.03, and i = 65 ° ±17° (presented in Fig. C.2). This optimization successfully fits the data across the NIRSpec IFU FoV (Fig. 2). Table 2 lists the best-fit parameters and associated errors.

Our adopting of a M/L_F814W gradient for $0\stackrel{″}{.}75<r<2\stackrel{″}{.}0$ is consistent with a color (F450W–F814W) change in the central region, which is likely caused by a change in the stellar population in the innermost regions, as is clearly visible in multicolor images of this galaxy (Fig. C3). This adoption of the M/L_F814W gradient stems from the fact that we did not calibrate the HST images for reddening in this work. This results in an improved fit to the variation in V_rms versus radius, and a tighter constraint on M_BH (see Figs. D2 and D3 for models that extend beyond the NIRSpec FoV). However, this constrained M/L_F814W is 3.5 and 13 times higher than what is predicted by the Salpeter initial mass function using PEGASE-2 models with solar metallicity and stellar ages of 13 Gyr (old) or 3 Gyr (intermediate), respectively; this resulted from the limited kinematic data and the absence of dust correction in our stellar mass model.

Figure E1 shows the 2D scatter plots for each parameter, with colored points representing their likelihood. While the other three parameters are well constrained by JAM, the inclination (i) is poorly constrained by the data (i > 47°) and is much higher than what was proposed from a morphological study (i = 35°; Erwin 2004). This could be explained by the slow rotation of the nucleus in NGC 4736, whose nearly spherical kinematics result in a similar appearance across a wide range of inclinations. Notably, the anisotropic orbital parameter constraints indicate that the nucleus of NGC 4736 is strongly dominated by tangential stellar orbits (β_r < 0).

An unusual feature of our JAM constraints is the positive covariance between M_BH and ${M/L}_{2\stackrel{″}{.}0}$ within the 3σ confidence level (CL) region of the 2D probability distribution function (PDF), where M_BH increases with ${M/L}_{2\stackrel{″}{.}0}$. This may be attributable to the continued rise in M/L_F814W beyond the NIRSpec FoV, which requires more extended kinematic data for accurate probing. In addition, the 2D PDF also reveals negative covariances between M_BH and ${M/L}_{2\stackrel{″}{.}0}$ with β_r at the same CL, indicating that the distribution of tangential stellar orbits becomes more prominent as M_BH and ${M/L}_{2\stackrel{″}{.}0}$ increase.

4. Discussion and summary

Given M_BH = 1.6 × 10⁷ M_⊙, a measured central velocity dispersion σ_c = 105 km s⁻¹, and a distance D = 4.66 ± 0.59 Mpc (Karachentsev et al. 2003), the SMBH’s SOI radius is calculated as R_SOI = GM_BH/σ_c² = 6.2 pc ($0\stackrel{″}{.}276$). This radius is two times larger than the FWHM of NIRSpec’s PSF, confirming that our stellar kinematic measurements are well resolved spatially. The kinematic signature of the SMBH is detected in the innermost ≈12 spaxels. Our M_BH measurement is the first in external galaxies⁶ accurately estimated using JWST observations and the stellar dynamical method.

We find that the model with M_BH = 2.5 × 10⁷ M_⊙ deviates from the JWST data by more than 1σ. We thus argue that M_BH = 1.6 × 10⁷ M_⊙ represents the most robust and reliable estimate. In particular, we have presented a model without a black hole along with models that assume M_BH = 10⁶ and 5 × 10⁶ M_⊙ (see Fig. 2). The model without a black hole and the one with M_BH = 10⁶ M_⊙ are nearly identical at the center and indistinguishable within the SOI. The model with M_BH = 5 × 10⁶ M_⊙ begins to rise above the lower-mass models, producing a distinct velocity signature of at least 5 km s⁻¹ at the center, which is more than three times larger than the measurement uncertainties (ΔV_rms). This difference is sufficient for the black hole’s presence to be detected and measured reliably.

We therefore conservatively conclude that M_BH,min = 5 × 10⁶ M_⊙ represents the lowest mass limit accessible with NIRSpec high-spectral-resolution observations using dynamical modeling. However, we note that this limit applies specifically to galaxies with central velocity dispersions σ_c ≲ 100 km s⁻¹ and distances D ≲ 5 Mpc.

Our M_BH constraint for NGC 4736 aligns well (within 1σ) with the M_BH–σ (Fig. 3; Kormendy & Ho 2013; van den Bosch 2016; Graham 2023b), M_BH–M_⋆,sph (Graham & Sahu 2023), and M_BH–M_⋆,gal (Graham 2024) correlations. While the relations from Kormendy & Ho (2013) and van den Bosch (2016) were derived by grouping spiral, lenticular, and elliptical galaxies, Graham (2023a,b) separated the galaxy types, and we used the spiral-galaxy-specific relation. We adopted a galaxy stellar mass of log(M_⋆,gal/M_⊙) = 10.51 ± 0.14 and a spherical-bulge stellar mass of log(M_⋆,sph/M_⊙) = 10.03 ± 0.21 (Graham & Sahu 2023) based on the galaxy-light decomposition and a diet-Salpeter initial mass function (Davis et al. 2019).

Fig. 3. Location of NGC 4736 with respect to the MBH–σ (left), MBH–M⋆,sph (middle), and MBH–M⋆,gal (right) scaling relations. Our new measurements are consistent within 1σ with the MBH–σ scaling relations for spiral galaxies.

Our inferred M_BH corresponds to a bolometric-to-Eddington luminosity ratio of L_bol/L_Edd ≈ 1.1 × 10⁻⁵, which we determined using the value of L_bol = 2.2 × 10⁴⁰ erg s⁻¹ derived by Constantin & Seth (2012) and rescaled to our adopted distance. This ratio places our M_BH within the range typical for LINER targets (Ho 2008). Our M_BH is consistent with the value derived from the updated FP of black hole activity (Gültekin et al. 2019), which is applicable to LLAGNs. This agreement resolves the longstanding debate on the M_BH estimates for NGC 4736, as our result aligns with the predictions from the M_BH–σ and M_BH–M_⋆ relations for spiral galaxies.

Acknowledgments

We thank the referee for useful comments and discussion that helped improve the quality of the manuscript. We would like to thank the ICISE staff for their enthusiastic support. This work is partially supported by a grant from the Simons Foundation to IFIRSE, ICISE (916424, N.H.). This research was funded by Vingroup Innovation Foundation (VINIF) under project code VINIF.2023.DA.057. RS acknowledges support from the INAF grant number 1.05.23.04.04. NT would like to acknowledge partial support from UKRI grant ST/X002322/1 for UK ELT Instrument Development at Oxford. IC’s research is supported by the SAO Telescope Data Center. He also acknowledges support from the NASA ADAP-22-0102 and HST-GO-16739 grants. MPS acknowledges support under grants RYC2021-033094-I, CNS2023-145506 and PID2023-146667NB-I00 funded by MCIN/AEI/10.13039/501100011033 and the European Union NextGenerationEU/PRTR. AWG thanks Bloomfield Higher Secondary School, Zunheboto, Nagaland, India, from where he worked. This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with program(s) #9042 and #9360. This work is based [in part] on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #02016.

References

Appendix A: Stellar spectral libraries

Stellar spectral libraries are essential tools for analyzing observed spectra and have broad applications in astrophysics. In this work, we employed two distinct types of stellar spectral libraries – one empirical, i.e., the X-shooter stellar library (XSL) and the other synthetic, i.e., the Extensive library of PHOENIX stellar atmospheres and synthetic spectra – to extract stellar kinematics from NIRSpec G235H/F170LP observations. As each type of library has its respective strengths and limitations, it is crucial to evaluate the robustness of our stellar kinematics extraction by comparing results derived from both.

Fig. A.1. Same as panels C, D, and E in Fig. 1, but showing kinematics derived with the PHOENIX high-spectral-resolution synthetic spectral template, providing a comparison to the results obtained with the XSL.

Empirical libraries provide spectra based on real observations, ensuring accurate representation of physical phenomena. For instance, absorption features due to atomic and molecular transitions are reproduced at the correct wavelengths with appropriate shapes and strengths. However, these libraries are inherently limited in spectral resolution and wavelength coverage. Extending empirical libraries often requires significant observational efforts and is constrained by instrumental capabilities. Moreover, the chemical compositions of the stars within such libraries are not precisely known, potentially introducing systematic uncertainties in abundance determinations. Additionally, the parameter space spanned by empirical libraries is often restricted, typically covering a limited range of effective temperatures, surface gravities, and metallicities.

In contrast, synthetic spectral libraries are constrained by the completeness of atomic and molecular line lists, uncertainties in line-broadening parameters, and assumptions in modeling, such as plane-parallel versus spherical geometries or local thermodynamic equilibrium (LTE) versus non-LTE conditions. Despite these limitations, synthetic libraries offer significant advantages: they can be tailored to span a wide range of stellar parameters, elemental abundances, wavelength regimes, and spectral resolutions. These flexible capabilities make synthetic libraries valuable for a variety of astrophysical analyses.

Our stellar kinematic derivations in the nucleus of NGC 4736, using two stellar spectral libraries (XSL and PHOENIX) with the NIRSpec G235H/F170LP IFU and PPXF (as detailed in Sect. 2), show consistency between the results within 5% when comparing the Light-of-Sight Velocity Distribution (LOSVD) maps in Fig. 1 (XSL) to those in Fig. A.1 (PHOENIX). However, it should be noted that for the subsequent analysis (i.e., SMBH mass estimates), we adopted the stellar kinematics derived using the XSL library as the canonical reference measurement, while those measured using the PHOENIX library were used for cross-checking.

We have demonstrated that our PPXF stellar kinematic measurements from the NIRSpec G235H/F170LP using the empirical stellar population template XSL, are consistent with those derived using the synthetic stellar population model PHOENIX. Previous studies have also reported similar consistency when deriving the nuclear stellar kinematics of the ultracompact dwarf galaxy M60-UCD1 (Seth et al. 2014) and a sample of nearby dwarf elliptical galaxies (Nguyen et al. 2018), based on Gemini/NIFS and VLT/SINFONI data, using the synthetic PHOENIX templates and the empirical Wallace & Hinkle (1996) library of high-resolution spectra of ordinary cool stars.

A.1. X-shooter stellar library (XSL)

We used the third data release (DR3⁷) of the X-shooter stellar library (Verro et al. 2022), which includes 830 stellar spectral templates from 683 stars. The DR3 spectra covers the full wavelength range of the X-shooter spectrograph, spanning 3000–25,000 Å, with a resolution of R ≈ 10, 000 and dust extinction corrected. The library includes a broad range of stellar types, from O to M, including asymptotic giant branch stars, and covers most of the Hertzsprung–Russell diagram.

A.2. Library of PHOENIX stellar synthetic spectra

We cross-checked our stellar kinematics measurements obtained with PPXF using an alternative high-resolution synthetic spectral library, PHOENIX⁸ (Husser et al. 2013) of 11 stellar templates. Although the templates cover a spectral range of 500–55,000 Å, when compared to XSL in the NIR, PHOENIX offers 50 times higher spectral resolution (R ≈ 500, 000), extends from 0.3 to 5 μm, and lacks gaps associated with telluric absorption bands. The library covers a parameter space of 2,300 K ≤T_eff ≤ 25, 000 K, 0.0 ≤ log(g)≤ + 6.0, −4.0 ≤ [Fe/H]≤ + 1.0, and −0.2 ≤ [α/Fe]≤ + 1.2. The synthetic stellar spectra were generated from the spherical mode of PHOENIX, which was employed to create model atmospheres, accounting for micro-turbulence in the stellar atmospheres.

Appendix B: Line spread function for NIRSpec

In Sect. 2 we adopt a constant Gaussian LSF with an instrument broadening of σ_LSF = velscale = 50 km s⁻¹ per pixel (Cappellari 2017), (or FWHM of 117 km s⁻¹ per pixel at 2.3 μm) to derive stellar kinematics via PPXF with the NIRSpec G235H/F170LP IFU (R = 2700). Here, we double check this assumption by re-deriving the LSF using the JWST/NIRSpec IFU spectrum of the standard star P330E (PID: 01538, PI: Gordon, Karl D.). P330E is a a solar-type star (G2V) with many absorption lines and has a relatively low velocity dispersion that’s necessary to measure the LSF.

The LSF may vary across the FoV of NIRSpec G235H F170LP, as has been observed in previous studies using Gemini/NIFS (Seth et al. 2010) and VLT/SINFONI (Nguyen et al. 2018). Here, we characterize potential variations in the LSF by fitting a high-resolution synthetic template (i.e., XSL or PHOENIX) to the JWST/NIRSpec spectrum of the calibration star P330E. This is done by convolving the template with a Gaussian kernel and performing a χ² minimization to determine the instrumental broadening. Given the context of deriving stellar kinematics using the CO band head region, we particularly fit one or more of the following absorption features: the ¹²CO(2–0) band head at 2.2935 μm, which is very strong in cooler stars; the Na I doublet at 2.2062 and 2.2090 μm, moderately strong in G–M stars and relatively sharp and isolated; the Ca I triplet at 2.261, 2.263, and 2.265 μm, commonly found in cool dwarfs and giants; and the Mg I line at 2.2814 μm, which can be blended but is sometimes strong depending on metallicity. Once they reach the spectrograph, these absorption features are dispersed such that their intensity distribution with wavelength reflects the instrument’s LSF. Assuming the central wavelengths of these absorption lines (corrected for small redshifts) are known, we treat their Gaussian profiles as representative of the LSF. After locating these absorptions, we define a region around each absorption, subtract the continuum, normalize the absorption flux to the deep value, and sum the profiles to characterize the LSF.

The spectral resolution across the detector has a median value of 8.52 Å FWHM at 2.3 μm (R = 2700), with values ranging from 7 to 11 Å FWHM across the G235H/F170LP wavelength range of 16,600–31,700 Å. The final step before extracting the stellar kinematics is to apply the same spatial binning to the LSF cube as was used for the final IFU of NGC 4736. This varying LSF is then incorporated into the kinematic extraction using PPXF. We found that using the spatially varying LSF, instead of the constant broadening (Sect. 2), did not significantly affect the results, as the velocity dispersion is relatively high.

Appendix C: Change in populations in the mass model

Fig. C.1. Comparison of the HST/WFPC2 F814W (left) and HST/NICMOS NIC3 F160W images (right) with their MGE models, shown in 2D surface brightness density. Black contours represent the data and red contours the model, highlighting their alignment at corresponding radii and contour levels.

Fig. C.2. Left: Comparison of the stellar kinematics measured from NIRSpec G235H/F170LP observations and the best-fit JAM (with 1σ error) with a constant M/LF814W (solid blue line) and M/LF160W (solid red line). Discrepancies between the data and model are apparent beyond $0\stackrel{″}{.}75$, as marked by the dash vertical line. Right: Best-fit radially varying M/LF814W profile (described in Sect. 2; solid blue line). It exhibits a step-like characteristic with $M/L_{0\stackrel{″}{.}75}=8.3$ M⊙/L⊙ at $r\le 0\stackrel{″}{.}75$ and $M/L_{2\stackrel{″}{.}0}=8.39{\mathrm{M}}_{\odot }$/L⊙ at $r\ge 2\stackrel{″}{.}0$. Between these radii, the M/LF814W profile increases linearly. Similarly, for the best-fit radially varying M/LF160W profile (Appendix D; solid red line), we adopt $M/L_{0\stackrel{″}{.}75}=3.38{\mathrm{M}}_{\odot }$/L⊙ at $r\le 0\stackrel{″}{.}75$ and $M/L_{2\stackrel{″}{.}0}=3.40{\mathrm{M}}_{\odot }$/L⊙ at $r\ge 2\stackrel{″}{.}0$. Here, the M/LF160W profile is vertically offset by a value of 4.95 M⊙/L⊙ for visualization.

Our dynamical model identified evidence for a subtle variation in the M/L_F814W profile within the annular region of 0$\stackrel{″}{.}$75<r<2$\stackrel{″}{.}$0, where the profile exhibits a linear increase (right panel of Fig. C.2). To investigate the origin of this variation, we examined the corresponding changes in the HST/WFPC2 (F450W–F814W) color. We constructed central 6″ × 6″ color map of NGC 4736, shown in Fig. C3, by converting the fluxes from the F450W and F814W images (measured in counts s⁻¹) into their respective 2D surface brightness maps (in units of mag arcsec⁻²).

To ensure the accuracy of the color map, we applied a cross-convolution technique, where each image was convolved with the PSF of the other filter (e.g., for the F450W–F814W color map, the F450W image was convolved with the F814W PSF, and vice versa). This approach mitigates artificial gradients near the galaxy’s center caused by differences in PSF widths between the filters. Subsequently, the background level for each image was determined in an annulus located at the maximum radius available in the observations (beyond 40″ from the nucleus) and subtracted from the images. Finally, the (F450W–F814W) color map was created using the ABmag-based zero points, with corrections applied for foreground extinction.

Variations in the color map clearly reveal signatures of two distinct and dominant stellar populations: younger stars in the nuclear region and older stars at larger radii. A transition in color is observed within the annular region of $0\stackrel{″}{.}75<r<2\stackrel{″}{.}0$, consistent with our finding of a linear increase in the M/L_F814W profile.

Furthermore, Fig. C3 shows that the color gradients extend clearly into the very central pixels (i.e., $r\lesssim 0\stackrel{″}{.}75$), indicating the presence of kinematically colder stellar populations toward the nucleus. This suggests that the appropriate radial M/L gradient should likely extend all the way to the center as well, which would, in turn, lead to a higher inferred M_BH. A more detailed investigation of this effect, which would require significantly deeper and more expensive observations, is beyond the scope of this paper and is planned for future work, using stellar population analysis based on Keck optical spectroscopy obtained with OSIRIS at a spatial scale of $\approx 0\stackrel{″}{.}035$ for the nucleus of NGC 4736.

Fig. C.3. Central 6″ × 6″ color map of NGC 4736, derived from HST/WFPC2 F450W and F814W images. We see a significant variation in the range $0\stackrel{″}{.}75<r<2\stackrel{″}{.}0$. This gradient is likely attributed to a transition in stellar populations within this region, supporting our adoption a radially varying M/LF814W profile in the galaxy mass model.

Appendix D: Mass model with the HST/F160W image

As we did not apply a reddening correction to the HST/WFPC2 F814W image used to derive the stellar-mass model in Sect. 2, this effect is likely accounted for the radially varying M/L_F814W profile defined in Sect. 3. We addressed this systematic effect by using HST imaging in multiple bands. In particular, we cross-checked our stellar-mass model using the HST/NICMOS NIC3 F160W broad-band image at 1.6 μm, which also has a pixel scale of 0$\stackrel{″}{.}$1. The data were taken on 2003 April 20 (Project ID 9360; PI: Kennicutt). They consist of four exposures of 96 seconds each, totaling 384 seconds.

We derived the stellar-light distribution of NGC 4736 from this NIR imaging data using the MGE technique, following the same procedure described in Sect. 2. In this process, we also modeled the PSF of the NICMOS NIC3 F160W filter using TINYTIM (Krist et al. 2011), then decomposed it with the image for the MGE. The parameters of each spatially deconvolved Gaussian component of the MGE model are listed in Table D.1, and the comparison with its F160W image is shown in the right panel of Fig. C.1. It is clear that there is less reddening in the F160W image than in the F814W image, as evidenced by the better agreement between the data and the MGE model across the FoV of the F160W image.

Figure D1 compares the cumulative luminosity, calculated from Table 1 for F814W and Table D.1 for F160W, and the corresponding stellar mass, derived from these luminosities using the constant M/L values M/L_F814W = 8.4 M_⊙/L_⊙ (from Sect. 3) and M/L_F160W = 3.4 M_⊙/L_⊙ (from Appendix E). The increased light in the NIR filter is compensated by the lower M/L returned from the dynamical models. However, there is more central mass within 2″ in the F160W mass model compared to that derived from the F814W mass model, resulting in a smaller M_BH.

As seen in the left panel of Fig. C.2, the constant M/L JAM models based on both stellar-mass models from the WFPC2 F814W and NICMOS F160W images fail to provide good fits beyond 0$\stackrel{″}{.}$75, although the former stellar-mass model shows a noticeable improvement in matching the data. We therefore applied the radially varying M/L_F160W profile described in Sect. 3, with $M/L(\sigma )={(M/L)}_{0\stackrel{″}{.}75}=3.38{\mathrm{M}}_{\odot }$/L_⊙ within $0\stackrel{″}{.}75$, where σ is taken from Table D.1. We show the best-fit M/L_F160W variation profile in the right panel of Fig. C.2, together with the same M/L_F814W profile, which is vertically offset by a constant value of 4.95 M_⊙/L_⊙ for visualization purposes. Although the linear increase of the M/L profile within the annulus $0\stackrel{″}{.}75<r<2\stackrel{″}{.}0$ is small, it does significantly improve the fit between the model and the data beyond $0\stackrel{″}{.}75$, as clearly seen by comparing the left panel of Fig. C.2 with Fig. D2. Moreover, the change in M_BH when switching the mass model from a constant M/L_F160W to one with a linear variation profile is negligible, indicating that reddening has a minimal effect on our dynamical modelling at NIR wavelengths.

In Fig. D2 we show the 1D V_rms profiles produced from the best-fit JAM models with varying M/L for both the F814W and F160W stellar mass models. Both models provide good fits to the stellar kinematic data across the NIRSpec FoV. However, the M_BH constrained from the F160W mass model is 19% lower than that derived from the F814W mass model, consistent within the 2σ uncertainty.

Figure D3 again presents the best-fit models with varying M/L, taken from Fig. D2, to demonstrate that they accurately reproduce the NIRSpec stellar kinematics data. We also show additional models without a black hole or with the minimum detectable black hole mass, M_BH,min = 5 × 10⁶ M_⊙ that could be accessible with NIRSpec high-spectral-resolution observations. This dynamical technique is strictly applicable to galaxies with central velocity dispersions σ_c ≲ 100 km s⁻¹ and distances D ≲ 5 Mpc.

Fig. D.1. Enclosed luminosity (left) and mass (right) profiles constructed from the stellar-light MGE models derived from the HST/WFPC2 F814W image, assuming M/LF814W = 8.4 M⊙/L⊙ (Table 1), and from the HST/NICMOS NIC3 F160W image, assuming M/LF160W = 3.4 M⊙/L⊙ (Table D.1). Profiles are plotted at their widest radial extends, i.e., ≈50″ for the WFPC2 F814W and ≈25″ for the NICMOS F160W images.

Fig. D.2. Comparison of the 1D Vrms profiles from the best-fit JAM model with MBH = 1.6 × 107 M⊙ and the radial M/LF814W variation profile (solid blue line), and an alternative model with MBH = 1.3 × 107 M⊙ and the radial M/LF160W variation profile (solid red line), within the radial extent of the NIRSpec FoV.

Fig. D.3. Comparison of the 1D Vrms profiles from the best-fit JAM model with radially varying M/L profiles and MBH = 1.6 × 107 M⊙ (using the F814W stellar-mass model; solid red line, $R_{\mathrm{SOI}}\approx 0\stackrel{″}{.}276$) or MBH = 1.3 × 107 M⊙ (using the F160W stellar-mass model; solid blue line, $R_{\mathrm{SOI}}\approx 0\stackrel{″}{.}226$) and an alternative model with MBH = 5 × 106 M⊙ (solid green line) against the JAM model with MBH = 0 M⊙ (solid black line) within a radial extent of 3″. Despite the small SOI of the SMBH ($R_{\mathrm{SOI}}\approx 0\stackrel{″}{.}226-0\stackrel{″}{.}276$), our best-fit MBH model exhibits a kinematic imprint extending across the NIRSpec FoV ($\approx 1\stackrel{″}{.}5$) compared to the JAM model without a black hole (which accounts only for the stellar mass component). This results in a systematic increase in Vrms by approximately ΔVrms ≈ 2 km s−1, with the shift becoming more pronounced for higher MBH values, as shown in Figs. 2 and D2. For reference, an SMBH at the lowest mass limit detectable by NIRSpec high-spectral-resolution observations using stellar dynamical techniques (MBH, min = 5 × 106 M⊙ claimed in this work; Sect. 4) does not produce a visible kinematic effect beyond $\approx 0\stackrel{″}{.}8$.

Appendix E: Jeans anisotropic model

The nucleus of NGC 4736 is identified as a typical spiral with kinematics best described using spherical projections. To model its dynamics, we compared the observed V_rms map with predictions from JAM with a spherically axisymmetry alignment of the velocity ellipsoid, based on the Jeans equations. This approach aims to estimate M_BH (sampled in logarithmic scale) and constrain additional linearly scaled parameters: orbital anisotropy (β_r), M/L_F814W (or M/L_F160W), and inclination (i). This JAM is well-suited for less rotating spirals (Cappellari 2008). It predicts the mean LOSVD, assuming an axially symmetric velocity ellipsoid aligned radially (σ_r ≠ σ_θ = σ_ϕ). This spherical configuration is specified by setting align=‘sph’ in the JAM_AXI_PROJ procedure of the JAMPY⁹ package (Cappellari 2020). The modelling accounts for the NIRSpec PSF and uses fixed parameter search ranges and initial guesses as detailed in Table 2.

We performed a Markov chain Monte Carlo (MCMC) simulation within the JAM framework to explore the best-fit parameter space (M_BH, M/L_{F814W or F160W} or ${M/L}_{2\stackrel{″}{.}0}$, β_r, i) and estimate the associated uncertainties. The adaptive Metropolis algorithm (Haario et al. 2001), implemented in a Bayesian framework using the ADAMET¹⁰ package (Cappellari et al. 2013), was utilized for this purpose. The MCMC chains comprised 3 × 10⁴ iterations, with the initial 20% discarded as a burn-in phase. The remaining 80% of the iterations were used to construct the full PDF. The best-fit parameters were identified as those corresponding to the highest likelihood within the PDF. Statistical uncertainties for all free parameters were calculated at the 1σ (16–84%) and 3σ (0.14–99.86%) CL (see Fig. E1).

For each stellar-light MGE model (i.e., either F814W in Table 1 or F160W in Table D.1), we conducted two separate MCMC optimization procedures, referred to as Optimization 1 and Optimization 2 (or Opt. 1 and Opt. 2 in Table 2 for the F814W MGE model and in Table E.1 for the F160W MGE model), for the JAM modelings. The first assumes a constant M/L_{F814W or F160W} profile, while the second incorporates a radially varying M/L_{F814W or F160W} profile with a linear gradient (Fig. E1 for its best-fit parameters constrained through MCMC).

The dynamical results presented in Sect. 3 for the F814W MGE model were derived using our adopted LOSVD obtained from the empirical XSL and the stellar mass listed in Table 1. To assess the potential impact of template choice on the SMBH mass measurement uncertainty for NGC 4736, we repeated a MCMC optimization using the LOSVD derived from the PHOENIX synthetic stellar spectra library (Appendix A.2 and Fig. A.1). This test yielded a totally consistent JAM model with M_BH =(1.55 ± 0.20)×10⁷ M_⊙, ${M/L}_{2\stackrel{″}{.}0}=8.41\pm 0.04$ M_⊙/L_⊙, orbital anisotropy β_r = −0.50 ± 0.03, and an inclination angle of i = 64 ° ±18° for the F814W mass model and M_BH =(1.28 ± 0.20)×10⁷ M_⊙, ${M/L}_{2\stackrel{″}{.}0}=3.41\pm 0.02$ M_⊙/L_⊙, orbital anisotropy β_r = −0.35 ± 0.03, and an inclination angle of i = 64.8 ° ±18.3° for the F160W mass model, suggesting the systemic related to the choice of stellar spectra libraries is minimal.

Fig. E.1. Posterior distributions for the radial M/LF814W (left) and M/LF160W (right) variation profiles obtained after the burn-in phase of the ADAMET MCMC optimization process for the JAM applied to the stellar kinematics of the nucleus of NGC 4736, measured using the NIRSpec G235H/F170LP observations, are presented. In each triangle plot, the analysis involves four parameters (i, MBH, $M/L_{2\stackrel{″}{.}0}$, and βr), shown as scatter plots for their 2D projected PDF with white indicates maximum likelihood and a CL at 1σ, while black corresponds to a confidence level at 3σ. The top histograms display the 1D projected distributions of the best-fitting model parameters, including the best-fit values (indicated by vertical black lines, corresponding to the highest likelihood within the PDF) and 3σ uncertainty (indicated by either side black dash lines). Insert plots in the top-right of each triangle plot illustrate the agreement by directly comparing the Vrms data with that generated from the best-fit JAM using the radially varying M/LF814W profile (left) and the constant M/LF160W value (right) in the same velocity range. It is evident that the constrained model shows excellent agreement with the data across the NIRSpec FoV.

All Tables

All Figures

Fig. 1. Stellar kinematic measurements extracted from the JWST/NIRSpec G235H/F170LP IFU. Panel A: Logarithmic integrated intensity collapsed along the spectral dimension (excluding the detector gap at ∼2.41–2.49 μm). Panel B: PPXF fit for the central VORONOI bin. Part of the central bin spectrum, which contains the stellar absorption features at ∼2.3 μm (CO band heads), is plotted in black, and the best-fit empirical X-shooter stellar library (XSL; Appendix A.1) is overlaid in red. The residuals (data-model) are illustrated in green. Two vertical gray lines indicate the spectral region where we actually fit the model templates to the spectra of all bins. Panels C, D, and E: Rotation (V), velocity dispersion (σ), and root-mean-squared velocity ($V_{\mathrm{rms}}=\sqrt{V^2+{\sigma }^2}$). White contours show the intensity, which decreases by 0.5 magnitudes from the center outward. We also performed the same measurements using the PHOENIX high-spectral-resolution synthetic spectral template (Appendix A.2), which are shown in Fig. A.1.

In the text

Fig. 2. Comparisons of the stellar kinematics measured from the NIRSpec G235H/F170LP observations and a range of JAMs. The black line shows the model with a constant M/LF814W, and colored lines the models that incorporate a radially varying M/LF814W plus 1σ error. The dashed vertical line marks the extent of the SMBH’s SOI. The data profile was obtained by averaging over circular annuli, and the associated errors are defined as the median of the Vrms uncertainties calculated by PPXF (i.e., ΔV and Δσ) of all bins within the circular annuli.

In the text

	Fig. 3. Location of NGC 4736 with respect to the MBH–σ (left), MBH–M⋆,sph (middle), and MBH–M⋆,gal (right) scaling relations. Our new measurements are consistent within 1σ with the MBH–σ scaling relations for spiral galaxies.
In the text

	Fig. A.1. Same as panels C, D, and E in Fig. 1, but showing kinematics derived with the PHOENIX high-spectral-resolution synthetic spectral template, providing a comparison to the results obtained with the XSL.
In the text

	Fig. C.1. Comparison of the HST/WFPC2 F814W (left) and HST/NICMOS NIC3 F160W images (right) with their MGE models, shown in 2D surface brightness density. Black contours represent the data and red contours the model, highlighting their alignment at corresponding radii and contour levels.
In the text

Fig. C.2. Left: Comparison of the stellar kinematics measured from NIRSpec G235H/F170LP observations and the best-fit JAM (with 1σ error) with a constant M/LF814W (solid blue line) and M/LF160W (solid red line). Discrepancies between the data and model are apparent beyond $0\stackrel{″}{.}75$, as marked by the dash vertical line. Right: Best-fit radially varying M/LF814W profile (described in Sect. 2; solid blue line). It exhibits a step-like characteristic with $M/L_{0\stackrel{″}{.}75}=8.3$ M⊙/L⊙ at $r\le 0\stackrel{″}{.}75$ and $M/L_{2\stackrel{″}{.}0}=8.39{\mathrm{M}}_{\odot }$/L⊙ at $r\ge 2\stackrel{″}{.}0$. Between these radii, the M/LF814W profile increases linearly. Similarly, for the best-fit radially varying M/LF160W profile (Appendix D; solid red line), we adopt $M/L_{0\stackrel{″}{.}75}=3.38{\mathrm{M}}_{\odot }$/L⊙ at $r\le 0\stackrel{″}{.}75$ and $M/L_{2\stackrel{″}{.}0}=3.40{\mathrm{M}}_{\odot }$/L⊙ at $r\ge 2\stackrel{″}{.}0$. Here, the M/LF160W profile is vertically offset by a value of 4.95 M⊙/L⊙ for visualization.

In the text

	Fig. C.3. Central 6″ × 6″ color map of NGC 4736, derived from HST/WFPC2 F450W and F814W images. We see a significant variation in the range $0\stackrel{″}{.}75<r<2\stackrel{″}{.}0$. This gradient is likely attributed to a transition in stellar populations within this region, supporting our adoption a radially varying M/LF814W profile in the galaxy mass model.
In the text

	Fig. D.1. Enclosed luminosity (left) and mass (right) profiles constructed from the stellar-light MGE models derived from the HST/WFPC2 F814W image, assuming M/LF814W = 8.4 M⊙/L⊙ (Table 1), and from the HST/NICMOS NIC3 F160W image, assuming M/LF160W = 3.4 M⊙/L⊙ (Table D.1). Profiles are plotted at their widest radial extends, i.e., ≈50″ for the WFPC2 F814W and ≈25″ for the NICMOS F160W images.
In the text

	Fig. D.2. Comparison of the 1D Vrms profiles from the best-fit JAM model with MBH = 1.6 × 107 M⊙ and the radial M/LF814W variation profile (solid blue line), and an alternative model with MBH = 1.3 × 107 M⊙ and the radial M/LF160W variation profile (solid red line), within the radial extent of the NIRSpec FoV.
In the text

Fig. D.3. Comparison of the 1D Vrms profiles from the best-fit JAM model with radially varying M/L profiles and MBH = 1.6 × 107 M⊙ (using the F814W stellar-mass model; solid red line, $R_{\mathrm{SOI}}\approx 0\stackrel{″}{.}276$) or MBH = 1.3 × 107 M⊙ (using the F160W stellar-mass model; solid blue line, $R_{\mathrm{SOI}}\approx 0\stackrel{″}{.}226$) and an alternative model with MBH = 5 × 106 M⊙ (solid green line) against the JAM model with MBH = 0 M⊙ (solid black line) within a radial extent of 3″. Despite the small SOI of the SMBH ($R_{\mathrm{SOI}}\approx 0\stackrel{″}{.}226-0\stackrel{″}{.}276$), our best-fit MBH model exhibits a kinematic imprint extending across the NIRSpec FoV ($\approx 1\stackrel{″}{.}5$) compared to the JAM model without a black hole (which accounts only for the stellar mass component). This results in a systematic increase in Vrms by approximately ΔVrms ≈ 2 km s−1, with the shift becoming more pronounced for higher MBH values, as shown in Figs. 2 and D2. For reference, an SMBH at the lowest mass limit detectable by NIRSpec high-spectral-resolution observations using stellar dynamical techniques (MBH, min = 5 × 106 M⊙ claimed in this work; Sect. 4) does not produce a visible kinematic effect beyond $\approx 0\stackrel{″}{.}8$.

In the text

Fig. E.1. Posterior distributions for the radial M/LF814W (left) and M/LF160W (right) variation profiles obtained after the burn-in phase of the ADAMET MCMC optimization process for the JAM applied to the stellar kinematics of the nucleus of NGC 4736, measured using the NIRSpec G235H/F170LP observations, are presented. In each triangle plot, the analysis involves four parameters (i, MBH, $M/L_{2\stackrel{″}{.}0}$, and βr), shown as scatter plots for their 2D projected PDF with white indicates maximum likelihood and a CL at 1σ, while black corresponds to a confidence level at 3σ. The top histograms display the 1D projected distributions of the best-fitting model parameters, including the best-fit values (indicated by vertical black lines, corresponding to the highest likelihood within the PDF) and 3σ uncertainty (indicated by either side black dash lines). Insert plots in the top-right of each triangle plot illustrate the agreement by directly comparing the Vrms data with that generated from the best-fit JAM using the radially varying M/LF814W profile (left) and the constant M/LF160W value (right) in the same velocity range. It is evident that the constrained model shows excellent agreement with the data across the NIRSpec FoV.

In the text