Improved Dynamical Constraints on the Mass of the Central Black Hole in NGC 404

The STIS spectroscopic observations (PID: 12611, PI: Seth) were taken on 2012 November 18 with the G430L grating and the $52\buildrel{\prime\prime}\over{.} 0$ × $0\buildrel{\prime\prime}\over{.} 1$ slit. This provides spectra over a wavelength range of 2900 Å to 5700 Å with a pixel size of 2.73 Å, and spectral resolving power of R ∼ 530–1040. Seven individual exposures were taken, each with an exposure time of ∼1890 s, for a total exposure time of 13230 s. The source was dithered along the slit and centered near the E1 aperture position to minimize charge-transfer inefficiency losses.

Reduced and rectified spectroscopic exposures (x2d files) were downloaded from the Hubble Legacy Archive (HLA). The STIS images suffer from numerous imaging defects, especially vertical defects extending over many rows. To remove these features from our dithered data, we created a median image without dithering which was subtracted from each individual, dithered image. Each dither was separated by nine pixels, thus, some background galaxy light will be included in this median image. However, this effect is very small; the maximum row-averaged flux of the median image is $\lt 25 \%$ of that at the outermost radius we analyze.

We then combined the seven individual exposures into a single two-dimensional (2D) spectrum. We applied integer offsets of nine pixels between each exposure after verifying these offsets using the position of the nucleus. Before combining exposures, we used sigma-clipping at the 3σ level to remove cosmic rays and bad pixels, and also flag bad pixels based on their data quality (DQ) values. Specifically, we included pixels with DQ $\lt$ 1000. We then combined the exposures by taking a mean value of all remaining unflagged and unclipped pixels from the seven exposure stack. We also created an error frame for the combined image by creating variance frames from the original error frames. Then for each pixel, we added the appropriate variance frames together and divided by the square of the total number of exposures.

The signal-to-noise (S/N) of the central pixel in the combined spectrum is 37 and 50 per pixel at 3700 and 5000 Å, respectively. The S/N drops off to 20 per pixel at 5000 Å for binned data at +(0475–0575) and −(0675–0825), which are the outermost data we use in our spectral analysis.

The HST imaging data we use in this work are taken from different cameras. They include WFC3 and WFPC2, ACS/WFC, and ACS/HRC imaging; the data we use in this work are summarized in Table 1. Reduced drizzled images were downloaded from the HST/HLA. The HST/WFC3 images were obtained contemporaneously with the above STIS spectroscopy observation, and we use these in combination with the spectroscopic data to create a 2D mass map. We discuss the creation of point-spread functions (PSFs) for the HST data in Section 2.4.

NGC 404 was observed with Gemini/NIFS on 2008 September 21–22 using the Altair laser guide star system. The nucleus was observed for a total of 4560 s (see S10 for details). These data were used to derive stellar kinematics (from the CO band-head region) and the molecular hydrogen kinematics (from the 2.12 μm H₂ 1-0S(1) emission line) by S10; we use these data in the stellar and gas dynamical modeling presented in this work in Section 4 and Section 5, respectively.

We briefly review the derivation of the stellar kinematics here; more details are found in S10 and Seth et al. (2008b). We use the penalized pixel fitting (pPXF) code of Cappellari & Emsellem (2004) to measure the stellar line of sight velocity distribution (LOSVD) including the radial velocity, velocity dispersion, and the non-Gaussian terms h₃ and h₄, which measure the skewness and kurtosis of the LOSVD (van der Marel & Franx 1993). We fit the strong CO bandheads in the wavelength range from 22850 to 23900 Å. Kinematics were derived using eight spectra from the templates of Wallace & Hinkle (1996); these include spectral types from G to M (which show significant CO absorption) and classes from the main sequence to supergiants. We measured the instrumental resolution of each spatial pixel using sky lines. The LOSVD errors were estimated via Monte Carlo simulations using a propagated noise spectrum. Errors in the radial velocities ranged from 0.5 to 6 km s⁻¹ changing with S/N (see Figure 8 of S10).

One feature of this data set that is important in our analysis is the assumed center of the NGC 404 NSC. S10 used both a dust-emission corrected photocenter at R.A. = 01^h09^m26999 and Decl. = 35°43${}^{\prime }$0491) for modeling the gas kinematics and a kinematic center at R.A. = 01^h09^m26999 and Decl. = 35°43${}^{\prime }$04$\buildrel{\prime\prime}\over{.} 89$ for modeling the stellar kinematics. The stellar kinematic center is about half a pixel (∼0025) south of the photocenter. The uncertainties on these measurements are ∼002. We will show in Section 5 that this uncertainty in center translates into large variations in our gas dynamical BH mass estimates, suggesting the gas dynamics results are not robust.

The different PSFs for our HST images at different filters and our Gemini/NIFS spectra is a complicating factor in this analysis. The HST PSFs play a significant role in creating the mass map of the nucleus while the Gemini/NIFS PSF is critical for the dynamical modeling. In this work, we use two different types of PSFs: (1) for the HST/WFC3 data, we use Tiny Tim PSFs (Krist 1995; Krist et al. 2011), and (2) for the Gemini/NIFS data we use a two-component inner Gaussian + outer Moffat profile PSF presented in S10.

For the HST/WFC3 PSFs, we create the model PSF for each WFC3 exposure using the Tiny Tim routine for each of the four individual flt exposures. The PSFs were created using the ${\mathtt{tiny}}{\mathtt{3}}$ task, which includes a charge diffusion kernel and the distortion of a PSF. We then insert these into mock flt images at the position of the nucleus in each individual exposure to simulate our observations. Finally, for each filter, we combine the four PSFs at each of the four dither positions into the final PSF using the Astrodrizzle package (Avila et al. 2012).

For the Gemini/NIFS images, we use the two component PSF fit to the NIFS data from S10. The fits are made by convolving an HST image of the NGC 404 nucleus to match the NIFS continuum image using a two component PSF: (1) an inner Gaussian component, and (2) a larger scale Moffat profile ${({\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}/{[(1+(r/{r}_{d})}^{2})]}^{4.765})$) whose size is constrained by fits to the telluric calibrators taken alongside the science images. S10 find the best fit model PSF with an inner Gaussian full width at half maximum (FWHM) of ∼012, and an outer Moffat with FWHM of ∼0 $\buildrel{\prime\prime}\over{.} 95$ profiles; each containing about half of the light.

Astrometry is an essential step in this work. Below we outline the astrometric alignments we carried out:

• 1.  
  Due to the small field of view of our WFC3 data, we obtain absolute astrometry for the images by aligning an ACS F814W image in which the NSC is saturated to the 2MASS point source catalog. The alignment uses 12 unsaturated sources with a root mean squared offset of $0\buildrel{\prime\prime}\over{.}$ 20; this represents our absolute astrometric accuracy. We then align our WFC3 F814W image to this image using bright compact stars around the core.
• 2.  
  Astrometric alignment of NIFS, WFPC2, ACS/HRC, and the new WFC3 data were then tied to the WFC3 F814W data. The centroid position of the nucleus in WFC3 F814W was used as the reference position to align each image. Although dust clearly affects the area around the nucleus, it appears that the center of the nucleus does not suffer from significant internal extinction as seen in the F547M–F814W color map (Figure 2 of S10) and the F336W–F814W color map (this work, Figure 5). The alignment on the nucleus can be verified by comparing the alignment of sources in these color maps with the alignment of the NIFS Brγ map and the HST Hα image (Figure 6 in S10).
• 3.  
  Finally, we align the STIS spectroscopy to the WFC3 image. To do this, we collapse the spectra to create a 1D image of the slit, multiplying by the appropriate filter curve to make synthetic F547M and F336W images. We then compare these 1D images to the astrometrically aligned WFC3 images in the same filters to obtain WCS information on the slit positions. Because of the low S/N in the STIS image at large radius, the fitting of the 1D images was performed at radii $\lt 0\buildrel{\prime\prime}\over{.} 8$. Both image comparisons gave consistent WCS coordinates.

In this section, we describe the procedure to characterize the possible AGN continuum contribution to the nuclear spectrum through examining the time variability of the nucleus in different bands. Nuclear variability has been suggested by several previous observations (Maoz et al. 2005; Seth et al. 2010), but here we provide the strongest evidence yet for broad-band nuclear variability in NGC 404. This variability is important both in providing additional evidence for an accreting BH in NGC 404, and for constraining our stellar M/L modeling.

The most compelling measurements to date have been the UV variability presented by Maoz et al. (2005). Comparing their 2002 HST/ACS HRC UV observations to 1994 HST/FOS spectra with an $0\buildrel{\prime\prime}\over{.} 86$ aperture (Maoz et al. 1998) and 1993 pre-correction HST/FOC imaging (Maoz et al. 2005), they suggest that the nuclear UV flux at 2500 Å in 2002 was a factor of ∼3 smaller than in the 1993–94 observations; this drop is consistent with the possible AGN continuum visible in their FOS spectra. Because that work had a number of similar data sets where fading was not seen, they argue that this drop in flux is robust despite the very different data sets and large apertures used in the comparison. Variable hot dust emission in the NIR was suggested by S10 based on the comparison of the 2008 NIFS K-band images with 1998 HST/NICMOS NIC2 observations in F160W. Within the central $0\buildrel{\prime\prime}\over{.} 2$, the NIFS K-band image is 30%–40% fainter than the NIC2 observations. This picture is consistent with the observed spectral signature of hot dust emission in the nucleus, but the comparison across bands and instruments leads to significant uncertainties (S10). The observations we present below provide more robust evidence for nuclear variability than any previous sets of observations.

We examine the variability of the nucleus using the new WFC3 images (all taken 2012 November 18) along with images in similar filters from WFPC2 and ACS/HRC (Table 1). There are repeat measurements in three filters: (1) our WFC3 F336W image can be compared to a previous ACS/HRC F330W image taken 2002 October 28 (Maoz et al. 2005), (2) our WFC3 F547M image can be compared to previous WFPC2 images in F555W (1996 October 21) and F555W (1995 September 06), and (3) our WFC3 F814W image can be compared to a previous WFPC2 F814W image (1995 September 06).

To eliminate mismatches in the background level, we first estimated the sky fluxes in an annulus of $10\buildrel{\prime\prime}\over{.} 0$–$13\buildrel{\prime\prime}\over{.} 0$ away from the nucleus (the maximum radius available in all observations). These sky levels were then subtracted off of all nuclear photometric measurements. Next, we performed aperture photometry in three radial ranges: the nucleus ($\lt 0\buildrel{\prime\prime}\over{.} 2$), and two annuli surrounding the nucleus, from $0\buildrel{\prime\prime}\over{.} 2$–$0\buildrel{\prime\prime}\over{.} 5$ and $0\buildrel{\prime\prime}\over{.} 5$–10. The photometry was converted into luminosities using appropriate zero-points for each observation.

We compute the ratio of luminosities between two epochs to gauge the level of nuclear variability. Figure 1 shows the nuclear luminosities ($\lt 0\buildrel{\prime\prime}\over{.} 2$) increase by 18%–25% in all three bands between 1994 and 2012. Specifically, nuclear luminosities increase ∼18% for F336W, ∼25% for F555W (F547M), and ∼24% for F814W band. For each filter, we choose one epoch to be the reference epoch (filled circles in Figure 1) and compare the other epochs to the reference epoch. We use a $0\buildrel{\prime\prime}\over{.} 2$ aperture to capture a majority of the nuclear flux; the encircled energy within $0\buildrel{\prime\prime}\over{.} 2$ is 90%, 87%, and 85% for the WFC3 F336W, F547M, and F814W filters, respectively. The outer two annuli (triangles and plus signs in Figure 1) can then be used to gauge the variations that might be caused due to bandpass and detector sensitivity. We find variations of $\lt 3 \%$ at all epochs for these two annuli; we note that formal uncertainties are much smaller than this. Thus, the nuclear variations we find are highly significant (6–8 times the largest variations seen in our annuli), providing the strongest evidence to date for broadband UV-NIR variability in the NGC 404 nucleus.

Zoom In Zoom Out Reset image size

The increase in nuclear UV (F336W) flux we see from 2002–2012 is opposite to the decrease in 2500 Å UV flux observed by Maoz et al. (2005) from 1993–2002; this argues against a single transient event (e.g., supernova) being responsible for the variability. Given that the increase in flux at all wavelengths is opposite to apparent decrease in the NIR observed by S10 over overlapping time periods, this suggests that the variability likely occurs on timescales smaller than the decade scales probed here. Finally, we note that the STIS data analyzed below were taken contemporaneously with the WFC3 data. The variability we see here is consistent with the ∼17% contribution of the AGN continuum at 3700 Å that provides the best-fit to the nuclear spectra (Section 3.2).

We fit a range of single stellar population (SSP) models to the nuclear STIS spectroscopic data to determine the ages of stars and M/L spatial variations within the NGC 404 NSC. We show that the nuclear spectrum is somewhat better fit with the addition of an AGN continuum component, and use this as our default model in the nuclear regions. The stellar population of the NSC shows evidence of a range of ages, but an intermediate age component (1–5 Gyr) is particularly dominant in the central $\lt 0\buildrel{\prime\prime}\over{.} 3$ of the NSC.

3.2.1. SSP Model Fitting Methodology

3.2.2. Nuclear Spectrum and Testing an AGN Continuum

Based on the evidence for variability we found in the previous section, it appears an accreting BH may contribute to the continuum flux near the center of the galaxy. We therefore test whether including a power-law AGN component provides a better fit to the STIS data of the nucleus. We fit the nuclear spectrum ($\leqslant 0\buildrel{\prime\prime}\over{.} 05$) using the SSPs described both with and without an AGN continuum component: ${f}_{\lambda }\sim {\lambda }^{-\alpha }$. To determine the best fit value of α, we fit the spectrum over the entire STIS wavelength range; α is varied from −4.0 to 4.0 in steps of 0.1, and the resulting is evaluated for each model. We find the best fit of $\alpha ={0.5}_{-0.3}^{+0.4};$ as shown in Figure 2, these values provide the best fit over the full wavelength range, and also in the UV. The left panel of Figure 2 shows the ${\chi }^{2}$ of the whole spectrum as well as the ${\chi }^{2}$ in just the UV part of the spectrum ($\lambda \lt 4000$ Å) where the AGN is expected to contribute the largest fraction of the flux.

Zoom In Zoom Out Reset image size

We then compare a fit with an AGN component with $\alpha =0.5$ to one without an AGN component and find that the inclusion of an AGN component provides a better fit. The residuals of these two fits are shown in the right panel of Figure 2; while both provide visually similar fits, some clear differences are seen. Most notably, the Balmer emission lines in the non-AGN case are stronger and show a similar strength at higher orders (an unphysical scenario). More quantitatively, the non-AGN Balmer line ratios are Hγ/${\rm{H}}\beta \sim 0.48$, Hδ/${\rm{H}}\beta \sim 0.34$, and ${\rm{H}}\epsilon$/${\rm{H}}\beta \sim 0.49$. In the AGN case, the higher order Balmer lines are weaker with Hγ/${\rm{H}}\beta \sim 0.47$, Hδ/${\rm{H}}\beta \sim 0.24$, and ${\rm{H}}\epsilon$/${\rm{H}}\beta \sim 0.15$. These ratios are consistent with Osterbrock (1989) (Case A recombination with T ∼ 5000 K) which gives Hγ/Hβ ∼ 0.458, Hδ/Hβ ∼ 0.250, and H$\epsilon$/Hβ ∼ 0.153. We note however that similar line-ratios to the non-AGN fits are seen at larger radii where a power law component would not be expected (e.g., Figure 3, right panel).

Zoom In Zoom Out Reset image size

We note that these emission lines were excluded from the fit and thus do not contribute to the ${\chi }^{2}$ value, but the left panels show the improved ${\chi }^{2}$ and residuals for the fit including an AGN component. To quantify how significant this difference is, we calculate the Akaike Information Criterion (Cavanaugh 1997; Burnham & Anderson 2002, p. 100) corrected for finite sample sizes (AICc). Given the best-fit ${\chi }^{2}$ value of 206 for the model with an AGN and ${\chi }^{2}$ of 209 for the model without an AGN, the AICc suggests that the AGN model is ∼5 times more probable than the no AGN model. While this alone is not compelling evidence for an AGN component, the somewhat better fits combined with the strong evidence for nuclear variability presented in the previous subsection lead us to include an AGN power-law component for our population synthesis fits at radii $\lt 0\buildrel{\prime\prime}\over{.} 2$.

Nuclear Line Emission: From the residual spectrum of the central three pixels ($\leqslant 0\buildrel{\prime\prime}\over{.} 05$ or 0.75 pc) of the STIS long-slit spectroscopy we fit the Balmer ${\rm{H}}\beta$ recombination emission line with a Gaussian and find a flux of $(2.8\pm 0.3)\,\times {10}^{-16}$ erg s⁻¹ cm⁻². In S10, the total nuclear Hβ flux in an aperture of $1\buildrel{\prime\prime}\over{.} 0\times 2\buildrel{\prime\prime}\over{.} 3$ was $\sim 1.08\times {10}^{-13}$ erg s⁻¹ cm⁻², thus the fraction of Hβ flux in the nuclear spectrum is only 0.2% of the total emission in the nuclear region.

From the nuclear ${\rm{H}}\beta$ flux we can estimate the nuclear ${\rm{H}}\alpha$ flux to use for comparison to the X-ray luminosity and to estimate a nuclear star formation rate (SFR). We assume the Balmer decrement measured by Osterbrock (1989) of Hα/Hβ ∼ 3.1 for the case of AGN, giving an Hα flux of $(8.72\pm 0.05)\times {10}^{-16}$ erg s⁻¹ cm⁻² and ${L}_{{\rm{H}}\alpha }=(9.9\pm 0.1)\,\times {10}^{35}$ erg s⁻¹ with Cardelli et al. (1989) interstellar extinction law. In the case of star formation, ${\rm{H}}\alpha$/Hβ ∼ 2.85, giving an Hα flux of $(7.9\pm 0.2)\times {10}^{-16}$ erg s⁻¹ cm⁻² and ${L}_{{\rm{H}}\alpha }=(7.8\pm 0.1)\times {10}^{35}$ erg s⁻¹.

Binder et al. (2011) detected an X-ray point source associated with the NSC with ${L}_{2-10\mathrm{keV}}={1.3}_{-0.5}^{+0.8}\,\times {10}^{37}$ erg s⁻¹. We compare the nuclear ${L}_{{\rm{H}}\alpha }$ to this X-ray luminosity to determine its consistency with AGN emission. From this Hα flux, we can predict the X-ray flux using the empirical relationship of Panessa et al. (2006), with revised fit by Nguyen et al. (2014). The observed X-ray luminosity is consistent with the predicted flux from this relationship (${L}_{2-10\mathrm{keV}}={8.7}_{-5.0}^{+58.4}\times {10}^{36}$ erg s⁻¹). The uncertainty in the prediction here is very large but shows that the emission line signal we see is consistent with an AGN source to the X-ray emission.

If we instead assume the emission is coming from star formation, we can use the ${\rm{H}}\alpha$ luminosity ${L}_{{\rm{H}}\alpha }=(7.8\pm 0.1)\,\times {10}^{35}$ erg s⁻¹ to estimate an upper limit on the nuclear SFR. Using the relation of Murphy et al. (2011), we find a nuclear SFR of $(4.2\pm 0.4)\times {10}^{-6}$ ${M}_{\odot }$ yr⁻¹. This SFR translates into an expected radio luminosity at 5 GHz of $(1.6\pm 0.2)\,\times {10}^{22}$ erg s⁻¹ Hz⁻¹ (Murphy et al. 2011); this is significantly lower than the observed total radio luminosity of $3.4\times {10}^{24}$ erg s⁻¹ Hz⁻¹ at the center of NGC 404 (K. Nyland et al. 2017, in preparation). The radio morphology shows a compact source consistent within astrometric errors of the photometric and kinematic center of our NIFS data; the unresolved fraction of the total flux is ∼75%, thus, this emission is still far too large to be related to the Hα emission via star formation. Thus, the nuclear Hα emission favors AGN-related nuclear X-ray and radio emission.

3.2.3. Stellar Population Modeling of STIS Data

To obtain sufficient S/N (>20) for stellar population modeling, we bin the STIS spectrum along the slit. Near the center, the individual pixels ($0\buildrel{\prime\prime}\over{.} 05$) exceed this S/N threshold, but the outermost bins (at $r=0\buildrel{\prime\prime}\over{.} 675-0\buildrel{\prime\prime}\over{.} 825$) extend over $0\buildrel{\prime\prime}\over{.} 15$. All of the bins analyzed here had S/N $\gt \,20$ at wavelengths longer than 5000 Å. Bins with lower S/N appeared to have significant detector artifacts making stellar population fits unreliable. In this subsection, we focus on our best-fit results—these include an AGN component in the central $0\buildrel{\prime\prime}\over{.} 2$ and use the BC03 models. We discuss the uncertainties in our modeling in the next section.

The best-fit model and residuals to the central pixel are illustrated in the left panel of Figure 3. The model is shown in red with the individual SSP components shown in blue. The details of these stellar populations are shown in Tables 2 and 3. The fit in the central pixel is dominated by intermediate age stars with 70% of the light in the 1 Gyr SSP, 5% of the light in the 2.5 and 5 Gyr SSP, 17% of the light in AGN (see Table 2). Only a small fraction of the light is in younger populations ($\lt 570\,\mathrm{Myr}$). The best-fit A_V is 0.45 (${A}_{V}=1.086\times {\tau }_{V}$) and the M/L values are 0.55, 0.36, and 0.29 in the V, I, and H bands, respectively. The right panel of this figure shows a lower S/N spectral bin spanning $+(0\buildrel{\prime\prime}\over{.} 275-0\buildrel{\prime\prime}\over{.} 375)$ with higher extinction; the best fit ${A}_{V}\sim 0.93$ and the M/L values are 0.67, 0.52, and 0.38 in the V, I, and H bands, respectively. The reduced ${\chi }^{2}$ is 0.98 for the central pixel and 2.23 for the outer bin; the decrease in the goodness-of-fit at larger radii is likely due to increased contributions from bad pixels. The details of SSP model fits to all the spectral bins are presented in Table 3.

Table 2.  Stellar Population Fits to the STIS Central Pixel

Model	Age (Myr)	AGN	1	10	50	100	286	570	1000	2500	5000	10,000	13,000
	Z		0.02	0.02	0.02	0.02	0.02	0.02	0.02	0.02	0.008	0.008	0.0001
AGN	l	0.170	⋯	⋯	⋯	⋯	0.026	0.029	0.634	0.152	⋯	⋯	⋯
	m	⋯	⋯	⋯	⋯	⋯	0.028	0.034	0.735	0.204	⋯	⋯	⋯
No AGN	l	⋯	⋯	0.031	0	0.108	⋯	⋯	0.743	0.118	⋯	⋯	⋯
	m	⋯	⋯	0.002	0	0.024	⋯	⋯	0.707	0.267	⋯	⋯	⋯

Note. Best-fit stellar population of the central bin of the STIS spectrum. The upper model includes an AGN continuum component with power-law slope $\alpha =0.5$, while the lower panel does not include an AGN component. The best-fit light and mass fractions are given for each component (AGN or BC03 SSP) in the model. Here, l (light fraction) and m (mass fraction) for each best-fit model.

Download table as:  ASCIITypeset image

Table 3.  Resolved Stellar Population Fitting Results

r	Frac.	AGN	1	50	100	286	570	1.0	2.5	5.0	${\tau }_{V}$	$M/{L}_{I}$	$M/{L}_{I}$	${\chi }_{{\rm{r}}}^{2}$	$M/{L}_{I}^{\star }$
$(^{\prime\prime} )$			(Myr)	(Myr)	(Myr)	(Myr)	(Myr)	(Gyr)	(Gyr)	(Gyr)	BC03	BC03	CB07	BC03	BC03
(1)	(2)	(3)	(4)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)	(16)	(17)
+(0.475–0.575)	l	×	⋯	⋯	⋯	0.034	0.046	0.920	⋯	⋯	0.98	0.85	0.89	1.71	⋯
	m	×	⋯	⋯	⋯	0.081	0.065	0.854	⋯	⋯	⋯	⋯	⋯	⋯	⋯
+(0.375–0.475)	l	×	0.030	0.025	⋯	0.283	0.058	0.602	⋯	⋯	0.86	0.52	0.51	1.25	⋯
	m	×	0.002	0.038	⋯	0.257	0.093	0.615	⋯	⋯	⋯	⋯	⋯	⋯	⋯
+(0.275–0.375)	l	×	0.052	⋯	0.073	⋯	0.108	0.715	0.152	0.185	0.75	1.10	1.06	1.19	⋯
	m	×	0.154	⋯	0.091	⋯	0.118	0.372	0.170	0.192	⋯	⋯	⋯	⋯	⋯
+(0.175–0.275)	l	×	0.050	⋯	0.111	⋯	⋯	0.745	⋯	0.092	0.63	0.57	0.53	1.05	⋯
	m	×	0.057	⋯	0.115	⋯	⋯	0.656	⋯	0.126	⋯	⋯	⋯	⋯	⋯
+(0.075–0.175)	l	0.136	⋯	⋯	⋯	0.028	0.054	0.668	0.023	0.043	0.55	0.53	0.55	1.04	0.57
	m	⋯	⋯	⋯	⋯	0.045	0.143	0.660	0.062	0.095	⋯	⋯	⋯	⋯	⋯
+(0.025–0.075)	l	0.140	⋯	⋯	⋯	0.032	0.070	0.659	0.101	0.005	0.51	0.50	0.53	1.03	0.54
	m	⋯	⋯	⋯	⋯	0.078	0.098	0.563	0.157	0.014	⋯	⋯	⋯	⋯	⋯
–0.025—+0.025	l	0.170	⋯	⋯	⋯	0.026	0.029	0.634	0.152	⋯	0.41	0.36	0.40	0.98	0.39
	m	⋯	⋯	⋯	⋯	0.028	0.034	0.735	0.203	⋯	⋯	⋯	⋯	⋯	⋯
–(0.025–0.075)	l	0.156	⋯	⋯	⋯	0.034	0.027	0.648	0.106	0.019	0.45	0.40	0.44	1.02	0.43
	m	⋯	⋯	⋯	⋯	0.030	0.025	0.738	0.202	0.005	⋯	⋯	⋯	⋯	⋯
–(0.075–0.175)	l	0.142	⋯	⋯	⋯	0.057	0.063	0.599	0.077	0.052	0.48	0.49	0.47	1.03	0.52
	m	⋯	⋯	⋯	⋯	0.068	0.094	0.686	0.125	0.083	⋯	⋯	⋯	⋯	⋯
–(0.175–0.275)	l	×	0.035	⋯	0.128	⋯	0.072	0.691	⋯	0.074	0.51	0.50	0.48	1.05	⋯
	m	×	0.021	⋯	0.105	⋯	0.093	0.619	⋯	0.162	⋯	⋯	⋯	⋯	⋯
–(0.275–0.375)	l	×	0.109	⋯	0.093	⋯	0.145	0.481	⋯	0.173	0.59	0.56	0.58	1.12	⋯
	m	×	0.005	⋯	0.174	⋯	0.297	0.386	⋯	0.368	⋯	⋯	⋯	⋯	⋯
–(0.375–0.475)	l	×	0.155	0.187	0.184	⋯	0.084	0.250	⋯	0.279	0.71	0.63	0.65	1.17	⋯
	m	×	0.094	0.092	0.073	⋯	0.148	0.261	⋯	0.332	⋯	⋯	⋯	⋯	⋯
–(0.475–0.575)	l	×	0.068	⋯	0.180	0.109	0.264	0.220	0.160	⋯	0.88	0.45	0.47	1.25	⋯
	m	×	0.009	⋯	0.171	0.076	0.267	0.375	0.102	⋯	⋯	⋯	⋯	⋯	⋯
–(0.575–0.675)	l	×	⋯	⋯	0.129	0.309	0.310	0.073	⋯	0.179	1.03	0.47	0.50	1.31	⋯
	m	×	⋯	⋯	0.010	0.315	0.386	0.129	⋯	0.259	⋯	⋯	⋯	⋯	⋯
–(0.675–0.825)	l	×	⋯	0.084	0.066	0.331	0.316	0.160	⋯	0.042	0.98	0.38	0.36	1.43	⋯
	m	×	⋯	0.110	0.021	0.163	0.250	0.161	⋯	0.286	⋯	⋯	⋯	⋯	⋯

Note. Best-fit stellar population models of binned STIS spectrum. Column 1: radii included in bin; particularly, the five central most data points are single pixels, while the bins at larger radii combine multiple pixels. We therefore use the notation of ±(r₁–r₂)'' to specify the size of those spectra bins. Column 2: l (light fraction) and m (mass fraction) for each best-fit model. Column 3: AGN component; ×means no AGN component was included in the fit. Columns 3–12: SSP model age light and mass fractions. Column 13: Attenuation of starlight by dust. Columns 14–15: M/L_I for each spectral fit using the BC03 and CB07 SSP models. Column 16: reduced ${\chi }^{2}$ of each bin associated with BC03 SSP model. Column 17: the best-fit M/L${}^{\star }$ assuming no AGN contribution for the five central spectral bins.

Download table as:  ASCIITypeset image

The left panel of Figure 4 shows the light-weighted and mass-weighted ages along the slit. The mean light-weighted age is ∼2.5 Gyr and the mean mass-weighted age is ∼4 Gyr. A large fraction of the light appears to be coming from a 1 Gyr population across the slit. At radii $\lesssim 0\buildrel{\prime\prime}\over{.} 3$ (≲6 pc) this stellar population contributes ∼70% of the light and mass; this drops somewhat at larger radii (the right panel of Figure 4), continuing the trend seen at larger radii (Bouchard et al. 2010). This mass fraction is higher than the mass fraction of ∼1 Gyr found in the NSC by S10 using the same template fitting method using ground-based optical spectroscopy (see also Cid Fernandes et al. 2004), suggesting that this population is concentrated near the center of the NSC. This central region corresponds to the "extra-light" counter-rotating component observed in S10 and provides additional evidence that this portion of the NSC was formed from externally accreted gas during a merger ∼1 Gyr ago. This suggests that gas from minor mergers can effectively accrete into the very center of the galaxy to form an NSC.

Zoom In Zoom Out Reset image size

3.2.4. Assessing Uncertainties in the Stellar Population Modeling

We examine the correlation between the WFC3 colors and STIS spectroscopic M/L, and use this to create a new color–M/L relation for the NGC 404 nucleus. This is a critical step in improving the BH mass estimate by more accurately modeling the mass distribution in the NGC 404 nucleus.

We create three WFC3 color maps of F336W–F547M, F336W–F814W, and F547M–F814W by taking the astrometrically aligned image pairs and cross-convolving each image with the PSF of the other (e.g., for the F336W–F547M color map, the F336W image was convoluted with the F547M PSF and vice versa). The cross-convolution is necessary to accurately assess color gradients, since different width PSFs can introduce spurious gradients near the center of galaxies. We then create color images using the Vega-based zero points and correcting for foreground extinction (${A}_{{\rm{F}}336{\rm{W}}}=0.318$, ${A}_{{\rm{F}}547{\rm{M}}}=0.194$, ${A}_{{\rm{F}}814{\rm{W}}}=0.114;$ Schlegel et al. 1998). Figure 5 shows the WFC3 F336W–F814W (approximatelyU–I) color map of the nucleus of NGC 404. This map shows redder regions resulting from dust extinction on the eastern side of the nucleus while the western half of the nucleus appears to have little internal extinction with bluer areas in this region consistent with younger ($\lesssim 500$ Myr) populations. This color map is consistent with the WFPC2 F547M–F814W color map in Figure 2 of S10 (which was not cross-convolved).

Zoom In Zoom Out Reset image size

Next, we match up the STIS slit to the WFC3 color images to enable comparison of the spectroscopic M/Ls with the broad-band colors. The images were aligned with the STIS slit as described in Section 2.5. The white lines in Figure 5 show the location of the STIS slit on the F336W–F814W color map. Over the region where spectroscopic M/L values have been determined, the F336W–F814W color varies by more than 1 mag due to dust and stellar population variations. To account for the spatial variation of M/L due to both stellar population and dust extinction, we use the effective M/L (M/L${}_{\mathrm{eff}}$) that includes dust extinction using the prescription of Charlot & Fall (2000) and the best-fit ${\tau }_{V}$ from our models. The I band M/L${}_{\mathrm{eff}}$ along the STIS slit is shown in the bottom panel of Figure 6, while the cross-convolved F336W–F814W color along the slit is shown in the top panel. The M/L${}_{\mathrm{eff}}$ is highest in the reddened regions on the eastern side of the nucleus. Both the color and M/L${}_{\mathrm{eff}}$ show a minimum at the center, likely due to the contribution of an AGN continuum component.

Zoom In Zoom Out Reset image size

Figure 7 shows the correlation of the WFC3 F336W–F814W color map with the spectroscopic M/L${}_{\mathrm{eff}}$ values. We fit this correlation to get a color–M/L relation; the best-fit is the red solid line. The idea behind this approach is the same as in Bell & de Jong (2001), Bell et al. (2003), and Zibetti et al. (2009), but appropriate to the exact populations present in the NGC 404 nucleus as determined using the STIS population synthesis fits. We fit a linear relation between the logarithm effective spectral M/L${}_{\mathrm{eff}}$ and the F336W–F814W color. The error in this relation is determined in two ways: first, we propagate the errors on the spectroscopic M/L determined above, and second, we determine bootstrap errors of the fit by refitting using random sampling with replacement. The first method yields significantly larger errors than our bootstrapped errors, which are minimal. We calculate the 1σ uncertainty on our color–M/L relation from our total error budget, and show these as the thin blue (steeper) and yellow (shallower) lines in Figure 7, respectively. We will examine the effect of our color–M/L relation uncertainties on the dynamical models in Section 4.

Zoom In Zoom Out Reset image size

Figure 7 also compares our best-fit color–M/L relation to the one predicted by Bell et al. (2003) for the Sloan i band based on the Sloan u − i color. Even accounting for filter transformations, the slope of the Bell et al. (2003) relation (purple thick solid thin line) is significantly flatter than our best-fit relation and predicts a much higher mass than our derived relation. We also plot a color–M/L relation calculated in our filters by Roediger & Courteau (2015, green line), also using BC03 models with a Chabrier IMF. This relation has a similar normalization (as expected given the similarities in the models used), but is significantly steeper than our relation. The differences in slope play the most critical role in affecting our dynamical models, as the overall normalization is scaled to fit the kinematics. The differences in slope of our relation from the Roediger & Courteau (2015) and the Bell et al. (2003) relation are likely due to the unusual SFH in the NGC 404 nucleus, and suggest that knowing the local SFH provides us with important information in mapping out the stellar mass distribution in the nucleus.

As we discussed in Section 3.2.2, the fraction of young stars in our fits near the center of the galaxy depends significantly on the allowed contribution of the AGN (∼17%) to the spectrum. To quantify the effect this has on our M/L, we re-fit the color–M/L relation using the non-AGN fits for the central spectral bins (thin black line in Figure 7). The M/L${}^{\star }{\rm{s}}$ from these fits are on average 7% larger than fits including the AGN (column 17 in Table 3), and this leads to a slightly shallower slope in the color–M/L relation, but the difference is within our 1σ uncertainties.

We describe the creation of a mass map and model for the NGC 404 nucleus in this section. The first step in this process is to produce an M/L${}_{\mathrm{eff}}$ map for the nucleus by applying an M/L${}_{\mathrm{eff}}$–color relation to the color map. The M/L${}_{\mathrm{eff}}$ map in F814W is shown in the left panel of Figure 8. This was created using the fitted correlation of the F814W M/L${}_{\mathrm{eff}}$ versus F336W–F814W.

Zoom In Zoom Out Reset image size

It is straightforward to obtain the mass map for the NGC 404 nucleus by multiplying the M/L${}_{\mathrm{eff}}$ map with the F814W luminosity map pixel by pixel. The right panel of Figure 8 shows the WFC3 F814W mass map of the nucleus of NGC 404. Because the mass map is a weighted version of the F814W luminosity map, its PSF should be well described by the F814W PSF. Comparison of the white (mass map) contours with the red (luminosity) contours shows that the mass distribution in the nucleus is more symmetric than the F814W light emission profile. This is because the M/L${}_{\mathrm{eff}}$ map is corrected for dust extinction on the eastern side of the nucleus. We note that the central M/L${}_{\mathrm{eff}}$ value from spectroscopy is very close to the best-fit relation, and thus we make no special correction to the map due to the presence of the AGN. The mass profile of the central regions of NGC 404 is shown in Figure 10.

To create a mass model that can be used in dynamical modeling, we need to deconvolve the effects of the PSF on the mass map. We do this using a galfit model with three Sérsic components, the inner two to fit the complicated NSC, and the outer to fit the galaxy bulge. Because our mass map extends out to only 25'', we constrain the shape of the outermost bulge Sérsic component to the values derived by S10: Sérsic index n₃ = 2.5, ${r}_{\mathrm{eff},3}$ = 45'' (675 pc), and position angle P.A. = 80°, while leaving the normalization and flattening of this component as free parameters. These parameters are then fit to the 2D mass map along with all parameters of the inner two Sérsic profiles using galfit. The best-fit parameters along with error estimates propagated from the M/L errors of the color–M/L relation are shown in Table 4. The mass map, mass model, and residuals are shown in Figure 9. The contours of mass density in each panel show both data (black) and model (white) at the same densities to highlight the regions of agreement and disagreement between data and model. These show that our mass model is more symmetric than our luminosity model. Our color–M/L relation has yielded a fairly symmetric mass distribution by up-weighting the M/L in the dusty regions east (left) of the nucleus, and down-weighting the M/L where there are young stellar populations to the west (right) of the nucleus. This increased symmetry shows the success of our color–M/L relation.

Zoom In Zoom Out Reset image size

Table 4.  Galfit Mass Map Model Profile in WFC3 F814W^a

Component	Label	Sérsic n	${r}_{\mathrm{eff}}$	${r}_{\mathrm{eff}}$	P.A.	b/a	M
			(pc)	$(^{\prime\prime} )$	(°)		(×106${M}_{\odot }$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
Central excess Sérsic	Central S	0.5 ± 0.1	1.6 ± 0.1	0.11 ± 0.01	20.7 ± 0.8	0.974 ± 0.006	3.4 ± 0.2
Inner Sérsic	Inner S	1.96 ± 0.17	20.1 ± 0.4	1.33 ± 0.04	60.0 ± 3.7	0.958 ± 0.004	10.1 ± 0.1
Outer Sérsic (Bulge)	Outer S	2.5	675	45	80	0.997 ± 0.003	844.2 ± 6.4

Notes. galfit mass map model parameters of WFC3 F814W mass map. Columns 1 and 2: component name and its label (see Figure 10). Column 3: component index profile. Columns 4 and 5: effective radius or half-light (half-mass in the case of mass map) of the component profile in pc and '', respectively. Column 6: position angle. Column 7: ratio of major axis and minor axis. Column 8: mass of each component.

^aParameter values without error bars mean that we fix those parameters during galfit.

Download table as:  ASCIITypeset image

Figure 10 shows the 1D radial mass profile of the mass map versus the galfit model with the three components shown individually; the residuals in the bottom panel are <10% over the inner region. This is similar to the level of residuals in the surface brightness fit presented in S10. Both model and data profiles were determined using the iraf task ellipse (Jedrzejewski 1987).

Zoom In Zoom Out Reset image size

To incorporate this mass model into our dynamical modeling, we need to create a multi-Gaussian expansion (MGE) of the model (Emsellem et al. 1994). We do this by decomposing the galfit model into individual MGEs using the mge_fit_1d code by Cappellari (2002); the result is shown in Table 5. We obtain very similar results by parametrizing the PSF using an MGE and then directly fitting the 2D mass map, but prefer fitting the galfit model because of the superior handling of the PSF (it allows one to input the actual PSF as a 2D image). We fit our MGE out to ∼25''.

Table 5.  MGE of the galfit Mass Density Map Model

j	I0	${\sigma }^{\prime }$	${q}^{\prime }$
	(${M}_{\odot }$/${\mathrm{pc}}^{2}$)	$(^{\prime\prime} )$
(1)	(2)	(3)	(4)
1	303071.	0.048	0.974
2	258200.	0.051	0.974
3	7525.6	0.073	0.974
4	4395.6	0.084	0.974
5	2643.0	0.121	0.974
6	2632.1	0.191	0.974
7	16183.5	0.251	0.974
8	4477.1	0.310	0.958
9	3432.3	0.940	0.958
10	1722.6	1.039	0.958
11	1353.21	1.372	0.958
12	948.75	5.265	0.997
13	506.594	9.919	0.997
14	204.936	25.019	0.997

Note. MGE fit for our best-fit model, created from a galfit model of the mass map shown in Figure 9. Column 1: component number. Column 2: central surface density (Gaussian amplitude). Column 3: Gaussian width (in arcseconds) along the major axis. Column 4: axial ratio.

Download table as:  ASCIITypeset image

We compare our mass MGE to the mass MGE created in S10 from the surface brightness profile without any M/L variations in Figure 11. The two mass profiles are consistent with each other to within 7% between $0\buildrel{\prime\prime}\over{.} 27$ and $25\buildrel{\prime\prime}\over{.} 00$ (4–375 pc). However, at the very center, the new model is about 10% lighter than the original best-fit mass model.

Zoom In Zoom Out Reset image size

To propagate the error in the color–M/L relation into our mass maps, we also create mass maps using the 1σ uncertainties on the color–M/L relation. We also created mass maps from other color–M/L relations for other filters (i.e., in F547M and F336W). In general, we find similar results from all maps except for those that we calculate an F336W M/L. In these cases, excess F336W light at the center is not fully compensated for by the color–M/L relations and excess mass is inferred due to the AGN (see column 7 of Table 6).We will discuss the effects these uncertainties have on our dynamical models in Section 4.3.

Table 6.  Jeans Anisotropic Modeling–Stellar Dynamical Modeling Results in Different Mass Map Filters and Colors

j	Color	Filter	Map	${M}_{\mathrm{BH}}$	${\beta }_{z}$	γ	i
				104 [${M}_{\odot }$]			[°]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
1	F336W–F814W	F336W	galfit	${4.0}_{-2.5}^{+3.5}$	${0.15}_{-0.10}^{+0.10}$	${1.895}_{-0.105}^{+0.045}$	${23.0}_{-2.0}^{+2.5}$
2		F547M	galfit	${6.5}_{-0.8}^{+1.5}$	${0.00}_{-0.10}^{+0.05}$	${0.950}_{-0.090}^{+0.055}$	${21.0}_{-1.0}^{+1.5}$
3		F814W	galfit	${7.0}_{-0.4}^{+1.7}$	${0.05}_{-0.15}^{+0.05}$	${0.890}_{-0.045}^{+0.060}$	${20.5}_{-2.0}^{+1.0}$
4	F336W–F547M	F336W	galfit	${3.7}_{-2.3}^{+3.1}$	${0.10}_{-0.10}^{+0.15}$	${1.970}_{-0.120}^{+0.075}$	${23.0}_{-2.5}^{+2.0}$
5		F547M	galfit	${7.2}_{-0.7}^{+1.6}$	$-{0.05}_{-0.10}^{+0.10}$	${0.935}_{-0.105}^{+0.090}$	${21.5}_{-1.5}^{+1.5}$
6		F814W	galfit	${6.8}_{-0.8}^{+1.7}$	${0.05}_{-0.05}^{+0.05}$	${0.905}_{-0.090}^{+0.075}$	${19.5}_{-2.0}^{+1.0}$
7	F547M–F814W	F336W	Direct	${4.3}_{-2.3}^{+3.4}$	${0.10}_{-0.15}^{+0.10}$	${1.910}_{-0.120}^{+0.105}$	${22.5}_{-1.5}^{+2.5}$
8		F547M	Direct	${7.5}_{-0.9}^{+2.0}$	$-{0.05}_{-0.05}^{+0.05}$	${0.920}_{-0.090}^{+0.075}$	${21.5}_{-2.0}^{+1.5}$
9		F814W	Direct	${7.2}_{-0.5}^{+1.4}$	${0.05}_{-0.10}^{+0.05}$	${0.860}_{-0.105}^{+0.045}$	${20.0}_{-1.5}^{+2.0}$

Note. Exploring the systemic uncertainty in our dynamical models using mass maps created with different filters and colors. Column 1: number of nuclear mass map. Column 2: color used to create the color–M/L relation. Column 3: HST filter used in creation of the mass map. Column 4: method used to create MGE; Galfit means that mass maps were fitted to a Galfit model which was then fit to an MGE model, "Direct" means that the MGE was created directly from the mass map. Columns 5, 6, 7, and 8: best-fit BH mass (${M}_{\mathrm{BH}}$), anisotropy (${\beta }_{z}$), scaling factor (γ), and inclination angle (i) of each best-fit model determined from Jeans models. All error bars are determined by propogating the uncertainty in the color–M/L relation through the dynamical models (e.g., the blue and yellow lines in Figure 7 and Figure 13).

Download table as:  ASCIITypeset image

We use the Jeans anisotropic modeling (JAM) method and IDL software⁸ as the dynamical model to calculate the mass of the central BH, ${M}_{\mathrm{BH}}$. The model relies on an axisymmetric solution of the Jeans equations incorporating orbital anisotropy (Cappellari 2008). The anisotropy is characterized by a single anisotropy parameter: ${\beta }_{z}=1-{\sigma }_{z}^{2}/{\sigma }_{R}^{2}$, where ${\sigma }_{R}$ and ${\sigma }_{z}$ are the velocity dispersion in the radial direction and z-direction in ellipsoid aligned cylindrical coordinates. The model calculates the predicted second velocity moment, ${V}_{\mathrm{rms}}=\sqrt{{V}^{2}+{\sigma }^{2}}$, where V is the mean stellar velocity and σ is the velocity dispersion based on an MGE model. We note we are using a mass model that differs from the luminosity profile of the galaxy; this difference is incorporated into the JAM model by using separate luminosity and mass profile MGEs; however both are axisymmetric and thus may not fully capture the variations in kinematics due to e.g., dusty regions. We also note that because of the proximity of NGC 404, the kinematic maps have some contribution from partially resolved stars that makes the dynamical maps appear less smooth, especially at larger radii. The model has four parameters: (1) the mass of a central BH, ${M}_{\mathrm{BH}}$, (2) the mass scaling factor γ; while for a normal JAM fit, this would be the dynamical M/L, in our case this parameter is the ratio between the dynamical M/L and the M/L predicted by our stellar population fitting (which assumes a Chabrier IMF), (3) the anisotropy, ${\beta }_{z}$, and (4) the inclination angle, i.

The models were run over a grid in these parameters, the ${M}_{\mathrm{BH}}$ range is 3 × 10⁴${M}_{\odot }$–3 × 10⁶ ${M}_{\odot }$ and is gridded in step of ${\rm{\Delta }}\mathrm{log}$ ${M}_{\mathrm{BH}}$ = 0.1, γ is gridded in steps of 0.015 between 0.50 and 2.25, ${\beta }_{z}$ ranges from −1 to +1 with a step of 0.05, and the inclination runs from the minimum value allowed by our MGEs, i = 17° up to 315 in step of 05.

We compute the models on a grid of ${\beta }_{z}$, ${M}_{\mathrm{BH}}$, γ, and i. At each grid point we evaluate the ${\chi }^{2}$ of the the predicted ${V}_{\mathrm{rms}}$ relative to the kinematic data. We run the dynamical modeling to fit these four parameters (${\beta }_{z}$, ${M}_{\mathrm{BH}}$, γ, i) in 2D using 920 kinematic data points (degrees of freedom, dof) within a radius of $\sim 1\buildrel{\prime\prime}\over{.} 5$ of the center. Figure 13 shows ${\chi }^{2}$ contours as a function of ${M}_{\mathrm{BH}}$ versus the anisotropy parameter ${\beta }_{z}$ (left), scaling factor γ (middle), and inclination angle i (right) after marginalizing over the other two parameters. The minimum reduced ${\chi }_{r}^{2}\sim 1.26$ is found at ${M}_{\mathrm{BH}}$ = $7.0\times {10}^{4}$ ${M}_{\odot }$, ${\beta }_{z}\sim 0.05$ (i.e., nearly isotropic), $\gamma =0.890$, and i = 20°. The solid contours show the $1\sigma$, $2\sigma$ (white), and $3\sigma$ (red) levels (or ${\rm{\Delta }}{\chi }^{2}=2.30$, 6.18, and 11.83) for two parameters; we choose the confidence limits for two parameters to accurately capture the uncertainties shown in our two parameter plots above after marginalizing over the third and the fourth parameter in each plot. The BH is only detected at the 1σ level. Given the restrictions JAM models place on the orbital freedom of the system (relative to e.g., Schwarzschild models), it is common to use 3σ limits in quoting BH masses (see Section 4.3.1 for more details). Therefore, we use the 3σ upper limit of ${M}_{\mathrm{BH}}$ $\lt \,1.5\times {10}^{5}$ ${M}_{\odot }$. This is similar to (but slightly larger than) the upper limit of the light-follows-mass model in S10 (using a ${\rm{\Delta }}{\chi }^{2}=9$ limit).

Figure 12 shows the 1D ${V}_{\mathrm{rms}}$ versus JAM predictions for our axisymmetric mass model with BHs with different masses. The ${V}_{\mathrm{rms}}$ values shown are bi-weight averaged over circular annuli. These data are compared with lines for different BHs; these show that BHs with ${M}_{\mathrm{BH}}$ ≳ 10⁵${M}_{\odot }$ do not provide a good fit to the data. Particularly, the innermost bin looks consistent (within the errors) with the 3 × 10⁵${M}_{\odot }$ model. Beyond ∼03, the data look consistent with anything $\lt 3$ × 10⁵${M}_{\odot }$, and beyond $\sim 0\buildrel{\prime\prime}\over{.} 7$ there is little difference between the models.

Zoom In Zoom Out Reset image size

The most significant change in comparison with S10 is the anisotropy parameter, ${\beta }_{z}$. Our best model is more isotropic (${\beta }_{z}={0.05}_{-0.15}^{+0.05}$) than the best fit of ${\beta }_{z}=0.5$ found in S10. Given the observed isotropy in other galaxy nuclei (e.g., Schödel et al. 2009; Seth et al. 2014), we interpret this as a sign of the success of our mass model (which incorporates variations in M/L${}_{\mathrm{eff}}$) correctly representing the mass distribution in this system.

This analysis does not include any gas mass in the nucleus; we show here that we expect this assumption to have a minimal effect. H i gas is present at large radii in NGC 404 but has a hole at the center (del Río et al. 2004). There is CO emission from the regions in and around the nucleus: assuming a distance of D = 3.06 Mpc, the total molecular gas mass is estimated to be 6 × 10⁶${M}_{\odot }$ (Wiklind & Henkel 1990). Most of this gas is to the NE of the nucleus, and only a few percent is within the central 10'', thus the amount within the region we are modeling is ≲10⁵${M}_{\odot }$. This gas mass is less than 1% of the stellar mass in the region we are modeling the kinematics, and based on the molecular hydrogen emission map (Section 5), this gas is likely spread widely across the nucleus; thus we do not expect this component to affect our M/L or BH mass estimates. We examine the kinematics from the warm molecular hydrogen emission in Section 5.

The confidence intervals in our analysis above are based on the kinematic measurement errors and do not include any systematic uncertainties in the mass model. In this section, we examine the mass model uncertainties by analyzing additional, independent mass model images.

The primary uncertainty in ${M}_{\mathrm{BH}}$ comes from the central stellar mass profile. We can get a rough uncertainty on the BH mass by examining the uncertainty in the total mass in the central pixel of our mass map. The central pixel in our mass map has a projected stellar mass of $5.1\times {10}^{5}$ ${M}_{\odot }$ (with an area of 1.8 pc²); the dominant uncertainty in this photometric estimate of the central mass is the uncertainty in the central M/L${}_{\mathrm{eff}}$. This uncertainty is 20% or ∼1 × 10⁵${M}_{\odot }$. We can expect our BH mass uncertainty to be comparable to the uncertainty in the mass in this central pixel. We now examine the systematic uncertainty due to our mass map using several methods. Overall, we find that the systematic irrors on our BH mass are all consistent with the 3σ upper limit of $1.5\times {10}^{5}$ ${M}_{\odot }$.

4.2.1. Color–M/L Relation Errors

We calculated the 1σ uncertainties on our color–M/L relation, shown as blue and yellow lines in Figure 7. We created mass maps and MGEs from these 1σ steeper and shallower color–M/L relations and ran a full grid of JAM models for both. Figure 13 show the resulting 3σ confidence levels. The effect on the BH mass from these models is minimal, probably because even with the errors, the color–M/L relations are quite similar at the bluer colors found near the nucleus. However, the larger changes in γ and i show that these parameters are quite sensitive to the assumed color–M/L relation, and that our formal errors do not capture the full extent of the uncertainty in these parameters.

Zoom In Zoom Out Reset image size

As noted, the 1σ errors in our relation still do not encompass the color–M/L relations derived in previous work. To test the effect of these more extreme color–M/L relations, we also derived the stellar mass map based on the Bell et al. (2003) relation (the purple thick solid line in Figure 7) to examine how much the slope in the color–M/L relation affects the constraints on the BH mass. The resulting best-fit BH mass is lower, ${M}_{\mathrm{BH}}$ ${}_{,\mathrm{Bell}03}$ = $3.3\times {10}^{4}$ ${M}_{\odot }$, but the most notable difference is the mass scaling factor of ${\gamma }_{\mathrm{Bell}03}=0.33;$ this factor of ∼3 is expected due to the higher normalization of the Bell et al. (2003) relation.

We also fit our source using the color–M/L relation derived without including an AGN component in the fit. As noted in Section 3.3, the central M/L${}^{\star }$ and NSC mass increase by 7%, which suggests we will derive a less massive BH at the center of NGC 404. Using this non-AGN color–M/L relation to create our mass model, we get JAM models with a best-fit BH mass of ${M}_{\mathrm{BH}}$ $=\,$ 6.3 × 10⁴${M}_{\odot }$ and $\gamma =0.95$ fixing the anisotropy to ${\beta }_{z}=0.05$ and inclination i = 20°. Therefore, the assumption of an AGN contribution has minimal effect on our results, with a best-fit BH mass 9% higher than the non-AGN case.

To summarize, the errors in our color–M/L relation suggest minimal changes in our best-fit BH mass, and all reasonable models are consistent with our 3σ BH mass upper limit of $1.5\times {10}^{5}$ ${M}_{\odot }$. However, substantial changes in the best-fit γ and inclination are seen using different color–M/L relations.

4.2.2. Mass Maps from Other Filters

4.2.3. IMF Variations

In this section, we discuss additional sources of uncertainties due to our dynamical modeling techniques and the assumed center.

4.3.1. Dynamical Modeling Uncertainties

4.3.2. Central Position Uncertainties

We model the H₂ 1-0 S(1) transition kinematics with a tilted ring model, similar to the modeling approach used in S10. In summary, our models assume that the H₂ emitting gas is rotating in thin rings on circular orbits around the center of the nucleus of NGC 404. At each radius, the velocity of the gas traces the enclosed mass in the spherical symmetry approximation (see e.g., Neumayer et al. 2007)

$$\begin{eqnarray}&&{v}_{c}^{2}(r)=\displaystyle \frac{{GM}(\lt r)}{r}.\end{eqnarray} \tag{ 1 }$$

Each ring of gas can have its own inclination and position angle. We constrain these parameters for rings in the outer parts by fitting ellipses to the velocity field with Kinemetry (Krajnović et al. 2006). Within $0\buildrel{\prime\prime}\over{.} 1$ of the center of the nucleus, where the PSF is more important than in the outer parts, we optimize for the best values. As was found by S10, there is a kinematic twist close to the center, where the position angle (P.A.) changes by 70^◦.

For our model we generate a spectral cube with the same spectral sampling as the NIFS IFU, but spatially oversampled by a factor of 10. For each ring, we distribute its flux over the cube according to its spatial distribution and the velocity along the ring to mimic our observations. We also spectrally convolve the flux of each ring by a Gaussian with a width of 20 km s⁻¹ to mimic the intrinsic dispersion of the disk. We then convolve each spectral slice with the oversampled NIFS PSF and resample the cube in order to compare it with our data.

Contrary to S10, we fit our dynamical models directly to the kinematic data without first symmetrizing them to enable more reliable data versus model comparisons. We calculate the ${\chi }^{2}$ for both the predicted velocity field and the predicted dispersion field. However, since the measured velocity field is more constraining than the dispersion field, we use only the ${\chi }^{2}$ of the velocity field to determine our best-fit model. To optimize our fit for the BH mass, we follow S10 in only using the central 12 × 12 pixels for our ${\chi }^{2}$ determination. Using a wider field for the calculation of ${\chi }^{2}$ gives qualitatively the same results.

S10 modeled the surface brightness of the H₂ gas as two exponential disks. In reality, the distribution of the gas is very irregular, as is visible in the H₂ intensity contours in Figure 14. We therefore deconvolve the line flux derived from the IFU data of the H₂ gas by fitting a number of Gaussians to the peaks in emission close to the center with galfit. This different flux model leads to slightly different gas dynamical models compared to those of S10. However, this improvement in the modeling did only marginally change the best-fit BH mass.

Zoom In Zoom Out Reset image size

S10 found a best-fit inclination of 37° for the gas disk. We allow the inclination to vary between 17° and 47°, in steps of 5°. As with the stellar population models, we include an additional mass scaling factor ${\rm{\Gamma }}/{{\rm{\Gamma }}}_{\mathrm{SP}}$, which was allowed to vary between 1.0 and 1.5 in steps of 0.05. We used BH masses between (0–$4.5)\times {10}^{5}$ ${M}_{\odot }$ in steps of $2.5\times {10}^{4}$ ${M}_{\odot }$.

The exact location of the center of NGC 404 is uncertain. The photometric center, which was derived in S10 by correcting the Gemini/NIFS continuum emission for the emission of hot dust, is offset by 0.37 pc ($0\buildrel{\prime\prime}\over{.} 025=0.5$ pixels) from the kinematic center. We therefore run our grid of models for two different centers.

Figure 15 shows the results of H₂ 1-0 S(1) gas kinematics modeling. When using the kinematic center, as was done by S10, we derive a BH with mass ∼${2.0}_{-0.75}^{+1.25}\times {10}^{5}$ ${M}_{\odot }$(3σ). This is two times lower than the gas dynamical mass estimate found in S10, but within the errors of that previous determination. However, when using the photometric center we find that we cannot dynamically confirm the presence of a BH, i.e., our models are consistent with no BH. The best-fit mass scaling factor of the cluster is 1.325 ± 0.075. This is significantly higher than the mass scaling factor from the stellar dynamical models (${0.890}_{-0.060}^{+0.045}$).

Zoom In Zoom Out Reset image size

The reduced ${\chi }_{r}^{2}$ of our fits are high: ${\chi }_{r}^{2}=109$ for the kinematic center, and ${\chi }_{r}^{2}=164$ for the photometric center. Although it is possible that we have underestimated the uncertainties on the velocity data, the most likely contribution to the high ${\chi }_{r}^{2}$ is systematic; the models simply do not provide a good fit to the data. In particular, large residuals in velocity are found close to the center of the cluster (Figure 14). It is not possible to remove this (positive) velocity peak by shifting the center. We note that this residual peak was also visible in the residual map of S10 (see their Figure 16).

Since shifting the center of our models by 0.4 pc (half a pixel) has a huge influence on the ${\chi }_{r}^{2}$ and can either introduce or take away the necessity of a BH in our models, we do not consider the 2 × 10⁵${M}_{\odot }$ for the mass of the BH a significant result, especially since the peak in the dispersion map of the H₂ gas does not coincide with either the kinematic or photometric center. Fitting the H₂ gas with two disks, one with position angle (P.A.) of 50° and one with a P.A. of $-20^\circ$ can reproduce qualitatively some of the features seen in both the velocity map (the double-lobed structure) and the off centered peak seen in the dispersion map. However, we find that it is still not possible to model the positive central peak in velocity. It is likely that part of the observed motion of the H₂ gas is therefore not tracing the potential, and might be either due to infalling or outflowing gas, or foreground gas at larger galactocentric radii. This may also explain the mismatch in the mass scaling factor between the best-fit gas dynamical and stellar dynamical models.

Previous work has already provided strong evidence for the presence of an accreting BH in the nucleus of NGC 404, and our observations of nuclear variability (Section 3.1) and blue continuum in the STIS spectrum (Section 3.2) have added to this evidence. Our dynamical constraints provide a robust upper limit of 1.5 × 10⁵${M}_{\odot }$ on this BH. We note that our upper limit makes the NGC 404 BH consistent with the radio non-detection of a compact core in NGC 404 by Paragi et al. (2014); they place a ∼3 × 10⁵${M}_{\odot }$ upper mass limit on the BH based on fundamental plane measurements.

NGC 404 is thus unique; it contains an MBH whose mass is dynamically constrained to be ≲10⁵${M}_{\odot }$. This alone can provide some evidence on the mechanism that formed BHs in the early universe (e.g., Volonteri 2010; Greene 2012). The dynamical upper limits on possible BHs in M33 and NGC 205 are an order of magnitude lower than for NGC 404 (Gebhardt et al. 2001; Merritt & Ferrarese 2001; Valluri et al. 2005), however, there is no evidence that BHs exist in either of these objects. The lowest mass BH with a dynamical measurement is in NGC 4395 (4 × 10⁵${M}_{\odot }$; den Brok et al. 2015). However, some similar BH estimates do exist for very faint galaxies detected as broad-line AGNs (Reines et al. 2013; Moran et al. 2014; Reines et al. 2014). The lowest mass BH with any kind of mass estimate was recently presented by Baldassare et al. (2015) who find a virial BH mass estimate from the broad ${\rm{H}}\alpha$ line of 5 × 10⁴${M}_{\odot }$.

BH mass detections have also been claimed in globular clusters, and these are typically between 10⁴⁻⁵${M}_{\odot }$(Gebhardt et al. 2005; Noyola et al. 2008; Lützgendorf et al. 2011). However, these detections remain controversial, and lower mass systems can be explained by clusters of dark remnants at the centers of clusters (Lützgendorf et al. 2013; den Brok et al. 2014).

With our best-fit BH mass and upper limit, we can examine NGC 404 in the context of galaxy scaling relations. For massive ETGs, there is a tight correlation between BH mass and photometric bulge mass estimates (e.g., Kormendy & Ho 2013; McConnell & Ma 2013) or dynamical bulge mass estimates (e.g., Häring & Rix 2004; Saglia et al. 2016). The bulge mass of NGC 404 was determined in S10; with ${M}_{H}=-19.62$ and an $M/{L}_{H}=0.63$ as determined from stellar population synthesis, the total photometric bulge mass for NGC 404 is $9.5\times {10}^{8}$ ${M}_{\odot }$. S10 find the Sérsic index of the bulge to be $n\sim 2.5$, and this high Sérsic index is commonly thought to indicate a "classical bulge" formed from merging (i.e., Figure 6 of Fisher & Drory 2010). The upper limit on ${M}_{\mathrm{BH}}$ for NGC 404 of $1.5\times {10}^{5}$ ${M}_{\odot }$ clearly lies below the bulge-mass–galaxy-mass relation seen in other, mostly massive, ETGs (see Figure 16). An increasing number of mostly late-type galaxies have been found to fall below the elliptical bulge-mass–BH-mass relation, including the precision maser masses (Greene et al. 2010; Läsker et al. 2016); these studies also emphasize the ambiguity of defining the bulge component in these systems and the large scatter seen in BH mass for similar bulge masses. A break in the nearly linear bulge mass relationship seen in massive ETGs to a steeper bulge mass dependence has been proposed for "Sérsic" galaxies (those without a central core) and spiral galaxies (Scott et al. 2013; Savorgnan et al. 2016). Our upper limit on NGC 404 is consistent with this proposed break in the relationship; the Scott et al. (2013) relation is shown along with other scaling relations in Figure 16. We also plot (1) all ETG BH mass estimates and ${M}_{\mathrm{Sph}.}\lt \,{10}^{10}$ ${M}_{\odot }$, and (2) all dynamical mass estimates for BHs < 10⁶ ${M}_{\odot }$.

Zoom In Zoom Out Reset image size

Assuming bulge dispersion velocity of 40 ± 3 km s⁻¹ (Barth et al. 2002), we find BH masses of 0.8 to 3 × 10⁵${M}_{\odot }$ using the Kormendy & Ho (2013) relation for all galaxies, largely consistent with our BH mass upper limit. We also note that the dispersion value measured by Barth et al. (2002) is from a $2\buildrel{\prime\prime}\over{.} 0$ × $3\buildrel{\prime\prime}\over{.} 7$ aperture; kinematic data presented by Bouchard et al. (2010) suggest a somewhat lower average dispersion (∼35 km s⁻¹) within the bulge effective radius; this would imply an even lower mass BH for NGC 404.

The relation between total galaxy mass and BH mass was recently explored by Reines & Volonteri (2015). While NGC 404 clearly has an extended disk (Williams et al. 2010), the total mass of the galaxy appears to be dominated by the bulge: the 2MASS LGA value for the ${M}_{H}=-19.70$ (Jarrett et al. 2003), M/L_H = 0.34 and A_H = 0.026 (S10), and thus the total mass of the galaxy is also ∼10⁹${M}_{\odot }$. Comparing this to Reines & Volonteri (2015), we find that NGC 404 falls toward the bottom edge of active dwarf galaxies with BH masses estimated from single-epoch measurements of a broad line AGN. These, in turn, are located well below the BH-mass–galaxy-mass relation defined by other quiescent BHs with dynamical mass determinations.

Recently, efforts have been made to incorporate the bulge radius or density in the ${M}_{\mathrm{BH}}$–galaxy scaling relations (Saglia et al. 2016; van den Bosch 2016). In particular, Saglia et al. (2016) have found that some of the scatter in the e.g., ${M}_{\mathrm{BH}}$–σ relation is correlated with the average bulge density and radius. NGC 404's bulge has an effective radius of 0.675 kpc, and a density within this radius of ${\rho }_{h}\sim 6\times {10}^{8}$ ${M}_{\odot }$ kpc⁻³ (Table 4). Using the 2D correlations for all ETGs and classical bulges (CorePowerEClass) from Saglia et al. (2016), the prediction from the ${M}_{\mathrm{BH}}$–σ–r_h relation is $3.7\times {10}^{5}$ ${M}_{\odot }$, while the prediction from the ${M}_{\mathrm{BH}}$–σ–${\rho }_{h}$ relation is $4.5\times {10}^{5}$ ${M}_{\odot }$. Our upper limit is lower than both these predictions by a factor of ∼3.

These findings fit in a scenario of co-evolution of BH and classical-bulge masses, where core ellipticals are the product of dry mergers of power-law bulges and power-law Es and bulges the result of (early) gas-rich mergers and of disk galaxies. In contrast, the (secular) growth of BHs is decoupled from the growth of their pseudo-bulge hosts, except when (gas) densities are high enough to trigger the feedback mechanism responsible for the existence of the correlations between massive BH and galaxy structural parameters.

With the best fitting BH mass value, ${M}_{\mathrm{BH}}$ = $7.0\times {10}^{4}$ ${M}_{\odot }$, its sphere of influence radius is calculated based on the dynamics: ${r}_{g}=G$ ${M}_{\mathrm{BH}}$ $/{\sigma }^{2}$, where σ is the stellar velocity dispersion of the bulge, G is the gravitational constant. We find ${r}_{g}=0.3\pm 0.1$ pc or $0\buildrel{\prime\prime}\over{.} 02\pm 0\buildrel{\prime\prime}\over{.} 01$. Alternatively, we also estimate the BH sphere of influence radius via our mass model, where r_g is the radius at which the enclosed stellar mass within that radius is equal to twice the BH mass (Merritt 2015). The latter definitions of r_g give its value of 0.35 ± 0.07 pc or $0\buildrel{\prime\prime}\over{.} 023\pm 0\buildrel{\prime\prime}\over{.} 018$. Both methods give a consistent value of the sphere of influence of the BH in NGC 404. The fact that these values are similar to the pixel size of Gemini/NIFS and WFC3, and smaller than their spatial resolution, further justifies the upper limit placed on the BH with our dynamical models.

If we assume the NSC is described by the inner two components of our galfit mass model, the stellar population mass estimate is ∼1.35 × 10⁷${M}_{\odot }$ (Table 4). From the JAM modeling, we find that the mass scaling factor is $\gamma ={0.890}_{-0.060}^{+0.045}$, thus, the dynamical mass is ∼(1.2 ± 0.2) × 10⁷${M}_{\odot }$, somewhat lower than this population mass estimate (which assumes a Chabrier IMF). This mass is roughly consistent with the $1.1\times {10}^{7}$ ${M}_{\odot }$ found for the NSC in S10; we note that a different definition of the NSC was used in S10—they did not fit the inner portion of the cluster, and the mass was calculated based on the best-fit M/L ratio multiplied by the luminosity of both the central un-fit component plus the fitted NSC component. The contrast is clear in comparing the effective radii; the S10 NSC has an effective radius of $0\buildrel{\prime\prime}\over{.} 74$ (11 pc) while the best-fit here has an effective radius of $\sim 1\buildrel{\prime\prime}\over{.} 0$ (15 pc), which is estimated from the central two Sérsic components combined (see Table 4). This latter value is on the large side, but within the distribution of NSC sizes (e.g., Georgiev & Böker 2014). Our NSC mass estimates are consistent with the prediction of ${M}_{\mathrm{NSC}}\mbox{--}{M}_{\mathrm{Gal},\mathrm{dyn}}$ scaling relationship of Scott & Graham (2013) and Ferrarese et al. (2006).

As was found in S10, the dynamical and stellar population mass estimates for the NGC 404 nucleus agree quite well, with a mass scaling factor $\gamma =0.890$. This provides some evidence that the IMF in the nucleus of NGC 404 is roughly consistent with the assumed Chabrier IMF. This relatively "light" IMF is consistent with the ratios of dynamical masses to stellar population model estimates in massive globular clusters in M31 (Strader et al. 2011), and in massive ultra-compact dwarf galaxies (UCDs) that are likely stripped nuclei (Mieske et al. 2013; Seth et al. 2014; C. Ahn et al. 2017, in preparation). There is growing evidence that the IMF is not universal among galaxies (e.g., Cappellari et al. 2012; Conroy & van Dokkum 2012; Cappellari et al. 2013; Kalirai et al. 2013; Spiniello et al. 2014); given the high stellar density of the NGC 404 nucleus, this suggests that a parameter other than SFR density may be driving these IMF variations.

Finally, we note again that the STIS data provide strong evidence for the dominance of a ∼1 Gyr old stellar population in the inner few parsecs of the NGC 404 NSC. Compared to the S10 stellar populations determined from a ground-based spectrum of the entire NSC, this 1 Gyr population is clearly enhanced in the center. S10 found that the cluster counter-rotates between this inner portion and the outer portion of the NSC; the new results suggest that the counter-rotating component may be a 1 Gyr old addition to the pre-existing NSC.