Six new supermassive black hole mass determinations from adaptive-optics assisted SINFONI observations

TL;DR

This study measures six new supermassive black hole masses in nearby early-type galaxies using adaptive optics and integral field spectroscopy, enhancing the understanding of black hole-galaxy co-evolution.

Contribution

It provides six new black hole mass measurements in intermediate-mass galaxies using advanced dynamical modeling techniques, expanding the current sample.

Findings

Black hole masses are consistent with MBH-σe scaling relations.

NGC 3640 and NGC 7049 have constrained lower black hole mass limits.

Considering stellar metallicity gradients reduces NGC 4570's black hole mass estimate by 30%.

Abstract

Different massive black hole mass - host galaxy scaling relations suggest that the growth of massive black holes is entangled with the evolution of their host galaxies. The number of measured black hole masses is still limited, and additional measurements are necessary to understand the underlying physics of this apparent co-evolution. We add six new black hole mass (MBH) measurements of nearby fast rotating early-type galaxies to the known black hole mass sample, namely NGC 584, NGC 2784, NGC 3640, NGC 4570, NGC 4281 and NGC 7049. Our target galaxies have effective velocity dispersions ({\sigma}e) between 170 and 245 km s^(-1), and thus this work provides additional insight into the black hole properties of intermediate-mass early-type galaxies. We combine high-resolution adaptive-optics SINFONI data with large-scale MUSE, VIMOS and SAURON data from ATLAS3D to derive two-dimensional stellar kinematics maps. We then build both Jeans Anisotropic Models and axisymmetric Schwarzschild models to measure the central black hole masses. Our Schwarzschild models provide black hole masses which are consistent with recent MBH-{\sigma}e scaling relations. NGC 3640 has a velocity dispersion dip and NGC 7049 a constant velocity dispersion in the center, but we can clearly constrain their lower black hole mass limit. We conclude our analysis with a test on NGC 4570 taking into account a variable mass-to-light ratio (M/L) when constructing dynamical models. When considering M/L variations linked mostly to radial changes in the stellar metallicity, we find that the dynamically determined black hole mass from NGC 4570 decreases by 30%. Further investigations are needed in the future to account for the impact of radial M/L gradients on dynamical modeling.

Tables

Table 1. Table 1: The sample

Table 2. Table 2: Details of the SINFONI observing runs

Galaxy	Date	PID	Pixel scale	N	N	$T_{exp}$	AO mode
			(mas)	of exp.	comb. exp.	(h)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
NGC 584	2007 Jul 23,24	079.B-0402(A)	100	3	3	0.75	NGS
NGC 2784	2007 Dec 12,29, 2008 Jan 01,02	080.B-0015(A)	100	9	9	4.5	NGS
NGC 3640	2010 Apr 08,09	085.B-0221(A)	100	13	12	3.16	LGS
	2013 May 08,09, May 09	085.B-0221(A)	100	8	7	3.16	LGS
	2013 Dec 27, 2014 Jan 07,24	291.B-5019(A)	100	14	12	3.16	LGS
NGC 4281	2010 Sep 04	085.B-0221(A)	100	4	3	2.25	LGS
	2012 Mar 20	085.B-0221(A)	100	4	4	2.25	LGS
	2013 May 09,11	091.B-0129(A)	100	19	16	2.25	LGS
NGC 4570	2013 May 07	091.B-0129(A)	100	16	15	2	LGS
NGC 7049	2005 Jun 08,09,10,14,19	075.B-0495(A)	100	20	16	6.25	NGS
	2005 Jun 27, Jul 02,03,07	075.B-0495(A)	100	12	9	6.25	NGS

Table 3. Table 3: SINFONI spatial resolution

Galaxy	FWHMN	FWHMB	fN	Strehl
	(arcsec)	(arcsec)
(1)	(2)	(3)	(4)	(5)
NGC 584	0.20 $\pm$ 0.02	0.74	0.54	13 %
NGC 2784	0.21 $\pm$ 0.02	0.50	0.74	11 %
NGC 3640	0.19 $\pm$ 0.02	0.56	0.41	14 %
NGC 4281	0.22 $\pm$ 0.04	0.90	0.86	10 %
NGC 4570	0.18 $\pm$ 0.02	0.58	0.47	15 %
NGC 7049	0.20 $\pm$ 0.03	0.61	0.67	13 %

Table 4. Table 4: HST archival data

Galaxy	PID	Instrument	Filter
(1)	(2)	(3)	(4)
NGC 584	6099	WFPC2	F555W
NGC 2784	8591	WFPC2	F547M
NGC 3640	6587	WFPC2	F555W
NGC 4281	5446	WFPC2	F606W
NGC 4570	6107	WFPC2	F555W
NGC 7049	9427	ACS	F814W

Table 5. Table 5: Summary of dynamical modelling results

	JAM					Schwarzschild
Galaxy	MBH	M/L	$\beta$	$i$	${\chi }^2/DOF$	MBH	M/L	${\chi }^2/DOF$	r${}_{\mathrm{SoI}}{}^{}/{\sigma }_{\mathrm{PSF}}$
	$(\times {10}^8$M${}_{☉}{}^{})$	(M☉/L${}_{☉}{}^{})$		(∘)		$(\times {10}^8$M${}_{☉}{}^{})$	(M☉/L${}_{☉}{}^{})$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
NGC 584	${1.93}^{\pm 0.06}$	${5.4}^{\pm 0.1}$	${0.05}^{\pm 0.01}$	${89}^{\pm 3}$	4.35	${1.34}^{\pm 0.49}$	5.4±0.2	0.99	1.6
NGC 2784	${1.69}^{\pm 0.1}$	${7.7}^{\pm 0.3}$	${0.04}^{\pm 0.04}$	77±13	1.86	${1.03}^{\pm 0.54}$	${6.7}^{\pm 0.7}$	1.08	1.8
NGC 3640	${0.99}^{\pm 0.13}$	${4.2}^{\pm 0.1}$	${0.14}^{\pm 0.04}$	85±9	13.90	${0.77}^{\pm 0.51}$	${4.2}^{\pm 0.2}$	1.65	1.1
NGC 4281	${4.91}^{\pm 0.15}$	${12.5}^{\pm 0.2}$	$-{0.03}^{\pm 0.02}$	75±9	9.93	${5.42}^{\pm 0.80}$	${9.3}^{\pm 0.3}$	3.98	2.1
NGC 4570	${1.19}^{\pm 0.09}$	${6.6}^{\pm 0.1}$	${0.22}^{\pm 0.02}$	74±1	3.98	${0.68}^{\pm 0.20}$	5.5±0.1	1.87	1.1
NGC 7049	${3.16}^{\pm 0.84}$	${11.4}^{\pm 0.4}$	${0.04}^{\pm 0.04}$	44±10	3.18	${3.17}^{\pm 0.84}$	${11.9}^{\pm 0.3}$	1.38	1.6

Table 6. Table 6: Schwarzschild dynamical modeling results with variable M/L

Galaxy	MBH [M${}_{☉}{}^{}]$	M/Ldyn/M/L∗	IMF
	$\left(M_{☉}\right)$
(1)	(2)	(3)	(4)
NGC 4570	${4.04}_{-1.2}^{+0.9}\times {10}^7$	${0.98}_{-0.02}^{+0.02}$	Salpeter
	${4.2}_{-0.4}^{+0.6}\times {10}^7$	${1.14}_{-0.03}^{+0.03}$	Kroupa rev

Table 7. Table 7: MGE parametrisation of the HST PSF

N584		N2784		N3640		N4281		N4570		N7049
WFPC2	F555W	WFPC2	F547M	WFPC2	F555W	WFPC2	F606W	WFPC2	F555W	ACS	F814W
norm	$\sigma$	norm	$\sigma$	norm	$\sigma$	norm	$\sigma$	norm	$\sigma$	norm	$\sigma$
	(arcsec)		(arcsec)		(arcsec)		(arcsec)		(arcsec)		(arcsec)
(1)	(2)	(1)	(2)	(1)	(2)	(1)	(2)	(1)	(2)	(1)	(2)
0.1972	0.0173	0.2430	0.0173	0.1951	0.0173	0.1995	0.0173	0.1956	0.0173	0.1517	0.0284
0.5832	0.0484	0.5591	0.0463	0.5847	0.0475	0.5631	0.0489	0.5863	0.0479	0.6483	0.0649
0.0923	0.1251	0.0902	0.1186	0.0960	0.1200	0.0498	0.1158	0.0927	0.1258	0.0983	0.1513
0.0676	0.3116	0.0685	0.3012	0.0687	0.3075	0.0680	0.1477	0.0668	0.3080	0.0620	0.4047
0.0186	0.4724	0.0391	0.8523	0.0144	0.5062	0.0782	0.3278	0.0172	0.4619	0.0166	0.8361
0.0409	0.8753	–	–	0.0411	0.8630	0.0414	0.8395	0.0414	0.8506	0.0231	1.6155

Table 8. Table 8: MGE parameters

	N584	(F555W)			N2784	(F547M)			N3640	(F555W)
j	log $M_j$	log $I_j$	${\sigma }_j$	$q_j$	log $M_j$	log $I_j$	${\sigma }_j$	$q_j$	log $M_j$	log $I_j$	${\sigma }_j$	$q_j$
	($M_{☉}$)	($L_{☉}p{c}^{-2}$)	(arcsec)		($M_{☉}$)	($L_{☉}p{c}^{-2}$)	(arcsec)		($M_{☉}$)	($L_{☉}p{c}^{-2}$)	(arcsec)
(1)	(2)	(3)	(4)	(5)	(2)	(3)	(4)	(5)	(2)	(3)	(4)	(5)
1	7.746	4.891	0.055	0.80	7.421	4.981	0.061	0.80	7.719	4.239	0.094	0.80
2	8.371	4.496	0.201	0.64	8.100	4.313	0.287	0.80	8.403	3.999	0.259	0.90
3	8.805	4.168	0.492	0.61	8.802	4.273	0.674	0.80	9.074	3.850	0.677	0.86
4	9.137	3.872	1.014	0.61	9.169	3.998	1.413	0.80	9.542	3.620	1.553	0.82
5	9.225	3.468	1.544	0.82	8.094	2.901	1.930	0.45	9.698	3.238	2.969	0.78
6	9.625	3.399	2.667	0.81	9.403	3.690	2.635	0.80	10.074	3.052	5.736	0.76
7	9.944	3.137	5.089	0.84	9.795	3.602	4.581	0.80	10.237	2.683	10.588	0.76
8	9.759	2.666	8.320	0.61	9.993	3.043	10.944	0.80	10.222	2.408	13.376	0.86
9	10.231	2.764	12.257	0.66	9.293	1.830	19.749	0.80	10.274	2.074	22.256	0.76
10	10.304	2.210	26.342	0.61	10.042	2.785	20.772	0.45	10.345	1.807	29.307	0.95
11	10.366	1.748	48.130	0.61	10.376	2.215	58.832	0.45	10.360	1.295	60.314	0.76
12	10.274	1.046	77.768	0.95	9.182	0.355	94.984	0.80	10.676	1.286	78.209	0.95
13	–	–	–	–	10.204	1.614	96.422	0.45	–	–	–	–

Equations

$$\log P(V_{rms}|i,M/L,{\beta }_z,\log M_{BH})\propto -\frac{1}{2}\sum _n\underset{{\chi }^2}{\underbrace{{\left(\frac{V_{rms}-⟨v_{los}^2{⟩}^{1/2}}{\delta V_{rms}}\right)}^2}}$$

Full text

1 Introduction

Most massive galaxies harbor a supermassive black hole (SMBH) in their centers. While black holes are invisible by their nature, their mass can be estimated using the motion of dynamical tracers (i.e., stars or gas) in combination with sophisticated dynamical models. The literature contains more than 100 robust dynamical black hole mass determinations, slowly growing into a statistically significant sample. Relating these measured black hole masses (MBH) to different host galaxy properties (such as bulge stellar mass, bulge velocity dispersion $\sigma_{e}$, Sérsic index n and star formation) revealed several noticeably tight correlations, i.e. MBH-L (Kormendy & Richstone 1995; Magorrian et al. 1998); MBH-$\sigma_{e}$ (Ferrarese & Merritt 2000; Gebhardt et al. 2000); MBH-n (Graham et al. 2001). Connecting vastly different scales these relations raise the question whether the growth of the black hole and the evolution of the host galaxy are entangled with each other (see recent reviews by Kormendy & Ho 2013 and Graham 2016). Current explanations suggest that black holes grow via two main processes: self-regulation by accretion of gas onto the black hole (facilitated by galaxy merging or accretion of gas) (Silk & Rees 1998; Fabian 1999; Di Matteo et al. 2008; Volonteri 2010) and by mergers of black holes (following dry major mergers). Kulier et al. (2015) and Yoo et al. (2007) show that accretion is the main channel of black hole growth, but galaxy mergers become relevant for more massive galaxies (see also Graham 2012; Graham & Scott 2013; Krajnović et al. 2018a). Based on the scaling relations we can see a clear trend that the more massive the galaxy is, the more massive is usually its central black hole. The exact shape of the various scaling relations is however still under debate. While early studies suggested a single-power law (Kormendy & Ho 2013), it is nowadays a question whether the fundamental relation between black hole and host galaxy properties scales as double power-law (Graham & Scott 2013) or has to be described by a three-parameter plane (van den Bosch 2016; Saglia et al. 2016). Moreover, Krajnović et al. (2018a) and Mezcua et al. (2018) recently reported an up-bending of the scaling relations with higher galaxy mass questioning the existence of one universal scaling relation. The search for a fundamental relation is made even more difficult by an increased internal scatter in both the low- and high-mass regime of the scaling relations. In order to understand and reduce the increased scatter, different observational strategies need to be developed. It is important to understand the different measurement methods with their associated systematic uncertainties by obtaining multiple MBH measurements with different methods for individual galaxies as was done, e.g. in Walsh et al. (2010); Barth et al. (2016); Davis et al. (2017, 2018c); Krajnović et al. (2018b). On the other hand, it is also important to figure out intrinsic scatter due to different galaxy formation scenarios by obtaining more and more homogeneous measurements over the complete SMBH mass range in order to strengthen current theories and ideas.

Our SMASHING sample (see Section 2 for details) was created to exploit the capabilities of natural and laser guide star adaptive optics (AO) systems at 8m ground-based telescopes. Its purpose is to fill up the scaling relations with additional MBH measurements of early-type galaxies. By the time of the creation of the project in 2009, the black hole mass measurements were almost exclusively populated by Hubble Space Telescope (HST) measurements with the exception of Nowak et al. (2008) and Krajnović et al. (2009) who pioneered a new method to measure MBH by using ground-based spectroscopy in combination with AO systems, laser and natural guide stars, respectively. The SMASHING survey was planned to expand the AO method to a large range of early-type galaxies with different velocity dispersions, from the low (100 km/s) to the high ($\approx$ 300 km/s) end. First results, based on observations with NIFS and GEMINI, were published in Thater et al. (2017) and Krajnović et al. (2018b). Unlike many other MBH measurements in the literature, we used both small (high spatial resolution) and large-field integral field spectroscopic (IFU) data for our measurements. The high-resolution kinematics are crucial to probe the orbital structure in the vicinity of the SMBH (also outside of its sphere of influence) and the large-scale kinematics are needed for constraining the global dynamical M/L, as well as to trace the influence of the stars on radial orbits, which pass close to the SMBH, but spend most of the time at large radii. Including both data sets provides more robust MBH measurements, especially if the sphere of influence is hardly resolved.

This paper is the first of a series of papers based on SINFONI observations. We present adaptive optics-assisted near-infrared integral-field spectroscopic data for six fast-rotating axisymmetric early-type galaxies to study the stellar kinematics in the vicinity of their central black hole. We begin by introducing the sample and its selection in Section 2, followed by the data acquisition and reduction in Section 3. In Section 4, we describe the extraction of the stellar kinematics from the near-infrared SINFONI and optical MUSE and VIMOS integral-field spectroscopic data. In addition to the kinematics, we combine high-resolution HST and Sloan Digital Sky Survey (SDSS) data to model the stellar surface brightness and thus examine the stellar brightness density of our target galaxies. In Section 5, we present the dynamical models which were constructed using two different and independent methods: Jeans Anisotropic Modeling (Cappellari 2008) and Schwarzschild (1979) orbit superposition modeling. We analyze our assumptions for the dynamical modeling with particular attention to M/L variations and discuss our results in the context of the MBH-$\sigma_{\mathrm{e}}$ relation in Section 6, and finally, conclude in Section 7.

2 The Sample

The six galaxies analyzed in this paper belong to our SMASHING galaxy sample to dynamically determine black hole masses in the nearby universe. Three of our target galaxies were selected from the ATLAS${}^{\textrm{3D}}$ volume-limited galaxy sample (Cappellari et al. 2011), from which one galaxy had already been observed in the SAURON project (de Zeeuw et al. 2002). The three remaining galaxies were observed with the VIMOS or MUSE instruments. Additional high spatial resolution data was obtained with the near-infrared SINFONI instrument to probe the direct vicinity of the SMBH. Based on their velocity dispersion, the sample galaxies are expected to be located in the intermediate $M_{\rm BH}$ range. The main properties of our six sample galaxies are summarized in Table 1.

Our target galaxies were selected based on a number of requirements for a successful $M_{\rm BH}$ determination. An important criterion for a robust black hole mass determination is the need to resolve the sphere-of-influence (SoI) of the black hole within which the SMBH dominates the galaxy potential. The SoI depends on the mass of the black hole $M_{BH}$ and the velocity dispersion of the galaxy within an effective radius $\sigma_{e}$ and is defined as $r_{\textrm{SoI}}=GM_{\rm BH}/\sigma_{e}^{2}$ where G is the gravitational constant. We calculated an estimated value for $r_{\textrm{SoI}}$ using black hole masses based on the $M_{\rm BH}-\sigma_{\rm e}$ relation from Tremaine et al. (2002)111The data acquisition process for this project started in 2008. In that time Tremaine et al. (2002) was one of the best representations of the black hole - host galaxy scaling relations. Tremaine et al. (2002) is very similar to the scaling relations that we show in Fig. 10. The selection based on the scaling relation by Tremaine et al. (2002) was only to select galaxies that were most likely to provide robust $M_{\mathrm{BH}}$ estimates. However, the required observing time and obtaining useful data in the near-infrared with LGS AO trimmed the sample more significantly than any scaling relation. and the ATLAS${}^{\textrm{3D}}$ velocity dispersions from Cappellari et al. (2013b). Using the large set of information from both large-scale and high-resolution IFUs, we can probe sphere-of-influences which are 2-3 times lower than the spatial resolution (Krajnović et al. 2009; Thater et al. 2017). With the goal to gain the best possible resolution, we utilized the AO mode from the SINFONI instrument using a natural guide star (if possible) or a laser guide star to correct for unstable seeing conditions.

Furthermore, archival HST imaging was needed for the galaxies of our sample to build detailed light models of the galaxy’s centers. We also ensured that the selected galaxies would not include any obvious bars or merger features indicating a non-relaxed galactic potential which would make the galaxies unsuitable for dynamical modeling with static potential models, as used here.

3 OBSERVATIONS

The mass measurement of massive black holes requires a large variety of data sets. Both high spatial resolution kinematic information of the central galaxy region to constrain the wide range of different stellar orbit families and large-scale IFU data to constrain the global galaxy characteristics are essential for a precise measurement. The IFU data is complemented by imaging data from HST and ground-based telescopes to construct a detailed mass model of the host galaxy. In the following section, we present the different observations from the IFUs towards the imaging data.

3.1 SINFONI IFU data

The near-infrared portion of our IFU data was obtained between 2005 and 2013 with the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI) instrument mounted on UT4 (Yepun) of the Very Large Telescope (VLT) at Cerro Paranal, Chile. SINFONI consists of the Spectrometer for Infrared Faint Field Imaging (SPIFFI) assisted by the adaptive optics (AO) module, Multi-Application Curvature Adaptive Optics (MACAO) (Eisenhauer et al. 2003; Bonnet et al. 2004). We observed each galaxy at K-band grating (1.94 - 2.45 $\mu$m) providing a spectral resolution of $R\sim 4000$ and a pixel scaling of 100 mas leading to a total field-of-view (FOV) on the sky of about $3.2\times 3.2\arcsec$ per pointing. Details of the observing runs for each galaxy are provided in Table 2. For each of our observations, we made use of the AO mode, either using a natural guide star (NGS) or an artificial sodium laser guide star (LGS) in order to correct for the ground-layer turbulence and achieve the highest spatial resolution possible. In the ideal case, the LGS mode still requires a natural guide star to correct for the tip-tilt disturbances in the wavefront, which are not sensed by the LGS. However, we often did not have a suitable tip-tilt star close to the galaxy and tip-tilt on the nucleus was not always possible, such that we applied the SINFONI Seeing Enhancer mode which provided a slight improvement to the natural seeing. Our observations show typical Strehl ratios of about 10 % (see Table 6). The observations were performed using the object-sky-object nodding scheme. At the beginning and end of each observing block, the respective standard star was observed at a similar airmass and with the same optical setup in order to correct for the telluric features at similar atmospheric and instrumental conditions. We used the SINFONI reduction pipeline to reduce the data and reconstruct the data cubes of the individual observations. This science frame contains spatial information in the X and Y directions and spectral information in the Z direction. As the data reduction was quite extensive, we mention a number of steps individually in the next subsections.

3.1.1 Data reduction and sky correction

The data reduction followed mostly the steps that are described in the SINFONI instrument handbook. The observations were reduced using the ESO SINFONI pipeline (version 2.4.8, Modigliani et al. 2007) in combination with additional external corrections to optimize the resulting data cubes. The ESO pipeline includes the bias-correction, dark-subtraction, flat-fielding, non-linearity correction, distortion correction and wavelength calibration (using a Neon arc lamp frame) for each observation of target and standard star. The nearest sky exposure was used to remove the night-sky OH airglow emission using the method described by Davies (2008). In the last step of the data reduction, each observation was reconstructed into a three-dimensional data cube.

3.1.2 Telluric and heliocentric velocity correction

A significant part of the data correction in the near-infrared regime is the correction for telluric absorption which originates in the Earth’s atmosphere (mainly ozone, gaseous oxygen, and water vapor). Telluric absorption lines are exceptionally deep at the blue end of the K-band and may vary over the time of the observation. Therefore, it is necessary to correct each science frame individually. Standard stars with known spectra are typically used to remove these atmospheric absorption features from science cubes.

For the telluric correction of the near-infrared spectra, we wrote a Python script to apply the same method as described in Krajnović et al. (2009). In most of the observation nights, two telluric stars were observed which gave us the opportunity to choose the telluric stars with a similar airmass to our science target. The telluric stars were either solar-like G0-2V stars or hotter B2-5V stars in an unsystematic order. We used the Python version 6.06 of the penalized Pixel fitting software333http://purl.org/cappellari/software (pPXF, Cappellari & Emsellem 2004) as upgraded in Cappellari (2017) to fit a stellar template showing the characteristic features of the telluric star. For the solar-like G-type stars, we use a high-resolution solar template (Livingston & Wallace 1991)444ftp://nsokp.nso.edu/pub/atlas/photatl/ and in the case of spectrally almost feature-less B-type stars we fitted a blackbody spectrum.

The telluric absorption corrected spectra were then corrected for the Doppler shifts due to the motion of the earth around the sun, commonly known as heliocentric correction. As some of our targets were observed at different times of the year, the velocity shifts in the observed spectra could be between 10-40 km/s, which is of the order of the LSF. We used our own python routine to correct the wavelength into the heliocentric frame of reference. The corrected wavelength is defined by $\lambda_{\textrm{corrected}}=(1+v_{\textrm{helio}}/c)\times\lambda_{\textrm{uncorrected}}$ where c is the speed of light and $v_{\textrm{helio}}$ is the projected heliocentric velocity which was calculated from the ESO pipeline for each data frame. The heliocentric correction was necessary for NGC 3640, NGC 4281 and NGC 7049, as in these cases the different observing blocks were spread widely throughout the year. We had to apply the heliocentric correction to each spaxel of each of our science frames individually.

3.1.3 Merging of the data cube

The individual frames of the observations from the different observing blocks were then combined spatially using the position of the center of the galaxy as reference. This position changed for each of our science frames as a dithering pattern of a few pixels was applied for each subsequent observation to ensure that the galaxy would not always fall onto the same pixels of the detector and thus adding systematic uncertainties. It was possible to identify the center of the galaxy with an accuracy of about one pixel (=0.05$\arcsec$) by comparing the isophotes of the reconstructed images and re-center them. In this step, we also excluded science frames with a bad PSF. Bad PSFs can be the result of poor seeing or an insufficient AO correction. In Table 2 we specify how many science frames were excluded for each galaxy. After the re-centering, we applied a sigma-clipping pixel reject algorithm to align the single science frames and created the final data cubes as in Krajnović et al. (2009) and Thater et al. (2017). The algorithm defines a new square pixel grid and interpolates the science frame to this grid. Flux values of the final data cubes were calculated as median flux values of the single data frames. Finally, we obtained $3\times 3\arcsec$ data cubes with 0.05$\arcsec$ pixel scale for the 100 mas SINFONI observations.

3.1.4 Correction of line-spread function inhomogeneities

In order to compare the spectra of the IFU with template spectra (which is needed for the extraction of the stellar kinematics, see Section 4.1), it is necessary to quantify the intrinsic dispersion of the SINFONI instrument. Therefore, we determined the spectral resolution of the SINFONI data from strong arc lines. While attempting to determine the spectral resolution of the SINFONI data, we encountered a problem: The spectral resolution over the full $64\times 64$ spaxel FoV was very inhomogeneous (see Figure 1) which was also recognized by Nguyen et al. (2018) and Voggel et al. (2018). In order to better characterize the shape and inhomogeneity of the LSF for the merged cubes, we applied the same data reduction routines to the respective arc lamp (except for the sky subtraction). We then built an arc line data cube by combining the reduced arc lamp frames. We used the same dither pattern as for the science frames to ensure that the arc line cube would fully resemble the data cube of the science object. From the combined arc line cube we then measured the LSF using six isolated, strong arc lamp lines for each spaxel. This LSF cube was later used when fitting each spaxel with a stellar template (details in Section 4.1). The spectral resolution across the field-of-view has a median value of $6.8$ Å  FWHM ($\lambda/\Delta\lambda$ = 3820) with values ranging from 5.5 Å  to 7.7 Å  FWHM.

3.1.5 Voronoi binning

The last step before determining the stellar kinematics was to ensure a sufficient and spatially uniform signal-to-noise (S/N) by spatially binning the final SINFONI data cubes with the adaptive Voronoi binning method555See footnote 2 (Cappellari & Copin 2003), Python version 3.1.0. In this method, based on the initial S/N estimate, spaxels are co-added while keeping the geometrical constraint of nearly round bins. An initially approximated noise estimate was obtained by median smoothing each spectrum with a kernel of 30 pixels width and calculating the standard deviation of the difference between the smoothed and the original spectrum. This initial estimate was then passed on as input S/N to the Voronoi binning script. The input S/N was systematically chosen between 50 and 70 balancing the desire to keep the central spaxels (if possible) unbinned to ensure a sufficiently high resolution in the center while increasing the quality of the outward spectra for the extraction of the kinematics. In our final binning scheme, we establish typical bin sizes of ¡ 0.1$\arcsec$ in the center, while 0.3-0.4$\arcsec$ diameter for bins at a radius larger than 1$\arcsec$.

3.1.6 SINFONI spatial resolution

The quality of our black hole mass measurements is indicated by the spatial resolution of the AO-corrected SINFONI data. We determined the spatial resolution by convolving high-resolution HST images with a double Gaussian model PSF and compared it to the collapsed image of the SINFONI IFU data. A detailed description of how we derived the spatial resolution is given in Appendix B.1. The resulting parameters are given in Table 6.

3.2 Large-field data

3.2.1 MUSE IFU data

The VLT/MUSE (Bacon et al. 2010) data of NGC 584 was taken on July 1st, 2016 under the science program 097.A-0366(B) (PI: Hamer). They obtained a total exposure time of 2700s divided into three $900$s on-source integrations each yielding a field of view of $60\times 60$ ″($\approx$ two effective radii of NGC 584). In addition, there was an off-source exposure of a blank field which can be used to estimate the sky. The fields were oriented such as to map the galaxy along the major axis with a large overlap, as every frame contained the nucleus of the galaxy. We reduced the data using the MUSE data reduction pipeline (Weilbacher et al. 2015), version 1.6. The reduction followed the standard steps, first producing master calibration files (bias, flat and skyflat), the trace tables, the wavelength solution and the line-spread function for each slice. Each on-target observation was reduced using these calibration files and closest in time illumination flats to account for temperature variations. In addition, a separate sky field and a standard star were reduced in the same way. From these, we extracted a sky spectrum and its continuum, as well as the flux response curve and an estimate of the telluric correction. The sky spectrum was applied to all three on-target frames, where we let the pipeline model the sky lines based on the input sky spectrum and the line-spread function. As all on-source frames contained the nucleus, we recorded its relative positions between the frames and applied the offsets with respect to the first one, prior to merging with MUSE pipeline merging procedure. In the final cube each pixel has the size of 0.2″$\times 0.2$″  and a spectral sampling of 1.25 Å  per pixel. For our MBH determination we only needed the high S/N central $30\arcsec\times 30\arcsec$ of the MUSE data cube and cut this region out. We then Voronoi-binned the cutted central region to a target S/N of 60 resulting into bin diameter sizes of $0.5\arcsec$ in the center and $3\arcsec$ at radii larger than 12″.

3.2.2 VIMOS IFU data

The large-field data for NGC 2784 and NGC 7049 were obtained between October 2006 and August 2007 using the VIsible Multi Object Spectrograph (VIMOS, Le Fèvre et al. 2003) mounted on UT3 Melipal under the science programs 078.B-0464(B) and 079.B-0402 (B) (PI: Cappellari).

The VIMOS data reduction was performed by Lagerholm et al. (2012) making use of the ESO pipeline777http://www.eso.org/sci/software/pipelines/ (version 2.3.3) and some IRAF tasks. It includes bias and sky subtraction, flatfield calibration, interpolation over bad pixels, cosmic-ray removal, spatial rectification, wavelength with HeArNe lamp exposures, flux calibration with standard stars and fringe-like pattern removal. As described in Lagerholm et al. (2012), they also corrected the fringe-like pattern in the spectral and the intensity variations in the imaging domain which were dominating the raw data. After the data reduction, they merged the individual science frames into final data cubes. In the same manner as for the SINFONI and MUSE data, we also Voronoi-binned the VIMOS data to a target S/N of 60, obtaining bin sizes of 0.5″in the galaxy center and 2-3″at radii larger than 7.5″.

3.2.3 SAURON IFU data

NGC 3640, NGC 4281 and NGC 4570 are part of the ATLAS${}^{\textrm{3D}}$ galaxy survey (Cappellari et al. 2011). The observations were obtained with the Spectrographic Areal Unit for Research on optical Nebulae IFU (SAURON, Bacon et al. 2001) at the 4.2-m William Herschel Telescope of the observatorio del Roque de los Muchachos on La Palma and reduced with the XSAURON software (Bacon et al. 2001) . The SAURON IFU has a FOV of $33\arcsec\times 41\arcsec$ with a sampling of $0.94\arcsec\times 0.94\arcsec$ square pixels, covering about 1-2 effective radii of our target galaxies. A detailed description of the stellar kinematics extraction of the ATLAS${}^{\textrm{3D}}$ sample is given in Cappellari et al. (2011). The resulting velocity maps of NGC 3640, NGC 4281 and NGC 4570 were already presented in Krajnović et al. (2011), and we show the full kinematic set of these galaxies in Figure 15. In addition, NGC 4570 is part of the SAURON survey (de Zeeuw et al. 2002) and presented in Emsellem et al. (2004). In this paper, we use the homogeneously reduced publicly available ATLAS${}^{\textrm{3D}}$ data888http://purl.org/atlas3d which was binned to a target S/N of 40.

3.3 Imaging data

For the high resolution central imaging of our target galaxies, we downloaded HST archival data. We obtained either Wide-Field Planetary Camera (WFPC2, Holtzman et al. 1995) or Advanced Camera for Survey (ACS, Ford et al. 1998) data from the ESA Hubble Science Archive, which generates automatically reduced and calibrated data. Cosmic rays were removed by taking the median of the aligned single CR-SPLIT images. Due to an unsuccessful sky subtraction in the archival data, the ACS image was reprocessed by applying the drizzlePac 999http://www.stsci.edu/hst/HST\_overview/drizzlepac package of the Astroconda distribution. For the large field of view imaging of our targets of the southern hemisphere, NGC 584, NGC 2784, NGC 7049, we used images of the Carnegie Irvine Galaxy Survey Project (Ho et al. 2011; Li et al. 2011; Huang et al. 2013). For the other three targets we used SDSS DR7 r-band images (Abazajian et al. 2009) which we received from the ATLAS${}^{\textrm{3D}}$ collaboration (Scott et al. 2013).

4 STELLAR KINEMATICS

4.1 Method

For each instrument, we independently measured the light-weighted stellar kinematics from the galaxy absorption line spectra using the Python implementation of the penalized Pixel Fitting method 111111See footnote 2 (pPXF, Cappellari & Emsellem 2004; Cappellari 2017). pPXF fits the galaxy spectrum by convolving a stellar spectrum template with the corresponding stellar line-of-sight velocity distribution (LOSVD), which is parametrized by Gauss-Hermite polynomials (Gerhard 1993; van der Marel & Franx 1993). In detail, the LOSVD is then specified by the mean velocity V, the velocity dispersion $\sigma$ and two additional quantities to describe asymmetric (h3) and symmetric (h4) deviations from a simple Gaussian. As the higher Gauss-Hermite polynomials are strongly coupled to the simple Gaussian moments, their relative weights are controlled by the so-called BIAS parameter which is dependent on the S/N of the data (Cappellari & Emsellem 2004; Emsellem et al. 2004). For low S/N data, the BIAS parameter prevents spurious solutions by biasing the derived LOSVD towards a simple Gaussian.

We analogously derived a second set of kinematics for each of our sample galaxies where we parametrized the LOSVD with the first two moments (V, $\sigma$) only. In this case, the BIAS keyword is not used by the code. This set of kinematics was needed to construct the Jeans Anisotropic Models (Section 5.2) which only take into account the lower-order moments of the LOSVD.

The usage of pPXF is twofold in order to minimize statistical variations across the field and reduce the computational expense. The first step is the creation of an optimal template by running pPXF on the global galaxy spectrum. The optimal template is a non-negative linear combination of the stellar library and consisted of typically 2-5 stars for the SINFONI data and about 30 stars for the large-scale data. Depending on the spectral range of the observed data, we used either MILES (Sánchez-Blázquez et al. 2006), Indo-US (Valdes et al. 2004) optical or Gemini Spectral Library of Near-IR Late-Type (Winge et al. 2009) stellar template library spectra, which are further described in the following two subsections. The optimal template is then used to fit the spectra from each Voronoi bin using $\chi^{2}$ minimization. While running pPXF on our spectra, we also added an additive polynomial to account for the underlying continuum. Furthermore, emission lines and regions of bad sky subtraction were masked during the fit. We then compared the fitted spectrum with the original spectrum for each bin. The standard deviation of the residuals (i.e., shown as green points in Fig. 2) was used to derive a final signal-to-residual-noise (S/rN) which measures both the quality of the data and the precision of the fit.

The errors of the recovered kinematics were derived with Monte Carlo simulations, the standard approach for LOSVD extractions. The complete measurement process is repeated for a large number of data realizations (500) where each realization is the original spectrum perturbed by adding noise in the order of the standard deviation of the pPXF residuals. Applying pPXF (with BIAS parameter set to zero) on each realization (with the same optimal template) provides 500 measurements of the LOSVD. The final error of each bin is then the standard deviation of the LOSVD distributions of these 500 realizations. The kinematic errors are spatially anti-correlated with the S/rN distribution, low in the center and larger in the outer regions. Mean errors are shown in Figure 4, where we compare the large and small-scale kinematics with each other. It is at first glance visible that the large-scale kinematics have much smaller errors than the SINFONI data (velocity: $\approx 2.5-5$ km s*-1* versus $\approx 5-10$ km s*-1* and velocity dispersion: $2.5-6$ km s*-1* versus $6-12$ km s*-1*).

4.2 SINFONI specifics

The SINFONI spectrograph in combination with adaptive optics provides spatially highly resolved spectra in the near-infrared regime yielding significant information about the motion of the stars surrounding the central black hole due its dust-transmissivity and its high resolution. A significant feature in the near-infrared is the CO absorption band head at about $2.3\,\mu$m which can be used to gain robust measurements of the LOSVD. We used the stellar template library by Winge et al. (2009) which consists of 23 late-type stars observed with the Gemini Near-Infrared Spectrograph (GNIRS) and 31 late-type stars observed with the Gemini Near-Infrared Integral Field Spectrometer (NIFS) to fit the SINFONI spectra in the range of 2.29 to $2.41\,\mu$m. Excluded from the fit were emission lines and incompletely reduced sky-lines which especially contaminated the third and fourth absorption line of the bandhead. Furthermore, to mitigate template mismatch effects in our kinematics extraction, we tested including both GNIRS and NIFS template stars as well as the restriction to only one instrument’s template stars. While all three attempts gave generally consistent results, the NIFS template stars could not always reproduce the Calcium line (at $\sim$ 2.25 $\mu$m) very well. This slight template mismatch often led to systematically lower velocity dispersions (but within the statistical errors). During the fitting procedure, we carefully examined and compared all three template library combinations and always chose the one that gave the best fit to the SINFONI spectra.

In order to recover reliable LOSVD measurements, we had to ensure that both the stellar templates and the SINFONI observations had a comparable spectral resolution before the fitting procedure. As the NIFS and GNIRS stellar spectra are provided at a better resolution than the SINFONI galaxy spectra, we had to degrade the template spectra to the same resolution as the SINFONI observations. Therefore, we convolved the template spectra ($\sigma_{\textrm{temp}}\approx 2.9-3.2$ Å) with a Gaussian having the dispersion of the difference between the dispersion of the Gaussian assumed LSF of the data and the stellar template. Our final pPXF fits reproduce the observed galaxy spectra very well as illustrated in Figure 2. For NGC 2784, NGC 3640 and NGC 4281 we also excluded the region around the Na I atomic absorption line at $\sim 2.2\,\mu$m as none of our stellar templates could match the line strength fully. As Silva et al. (2008) point out this is an often seen discrepancy between pure old galaxies and Galactic open cluster stars. We extended the masked regions because the blue part of the spectrum is very noise-polluted and biases the kinematics to a more noisy solution. Including or excluding this region, the changes in the four moments stay within the derived kinematical errors.

4.3 VIMOS & MUSE specifics

The kinematic extraction of the optical VIMOS and MUSE data was performed similarly to the ATLAS3D kinematic extraction. We re-extracted the kinematics for VIMOS data as the Lagerholm et al. (2012) extraction did not contain kinematic errors for each spaxel.

The optical IFU data matching stellar templates were taken from the medium-resolution Isaac Newton Telescope library of empirical spectra (MILES, Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011) stellar library121212http://miles.iac.es/ (version 9.1). We used the full sample consisting of 980 stars that span the wavelength range 4760-7400 Å  and fitted the wavelength range from 3800-6500 Å  in the galaxy spectra. As already mentioned, we also had to ensure that the instrumental resolution of the stellar templates and the optical data had comparable values. Beifiori et al. (2011) and Falcón-Barroso et al. (2011) report the instrumental dispersion of the MILES template library to be $\sigma_{\textrm{MILES}}=2.51$ Å. The MUSE spectral resolution was carefully measured by Guérou et al. (2017) based on sky emission lines, and the authors found a variation of the LSF with wavelength. In the for our analysis relevant wavelength range of 4800 to 6800 $\AA$ the spectral resolution changes from 2.5 to 2.9 $\AA$. In order to test the significance of this spectral resolution variation on our stellar kinematics measurement, we did an extraction using the extreme values of 2.5 $\AA$ and 2.9 $\AA$. The velocity dispersion changed in average by only 5 km/s which is within the kinematic error range, and we decided not to downgrade the MILES template library to the MUSE resolution. On the other hand, the VIMOS data have an instrumental dispersion of $\sigma_{\textrm{VIMOS}}=2.1$ Å (Rawle et al. 2008). Theoretically, a downgrading of the observed galaxy spectrum would be necessary here. As both instrumental spectral resolutions were relatively similar and we did not want to smooth relevant kinematic information, we did not convolve the VIMOS spectra to the lower resolution. All velocity dispersions are $\sigma>150$ km/s, so the effects of the slightly different resolutions on the derived kinematics are negligible. However, just like in Thater et al. (2017) we tested the effect of not downgrading the VIMOS spectra by using a well-defined sub-sample of the Indo-US stellar library (Valdes et al. 2004) as stellar template in the pPXF fit. Our sub-sample of the Indo-US stellar library consists of 52 spectra and covers a wavelength range of 3460 to 9464 Å  at a spectral resolution of $\sigma_{\textrm{Indo-US}}=1.35$ Å  (Beifiori et al. 2011). A comparison between the Indo-US kinematic extractions with the MILES extractions showed that the extracted kinematical maps displayed the same general features and trends. We could, however, discern a difference in the extracted values between the two stellar templates with the MILES velocity dispersions being systematically lower (10-20 km/s) and thus, more consistent with the stellar kinematics extraction from SINFONI. We furthermore recognized that the spectral fits were worse for the Indo-US fits such that we expected a template mismatch from the relatively small number of Indo-US template stars. Comparing our kinematic extraction with the extraction by Lagerholm et al. (2012) proved consistent results. We, therefore, decided to keep the MILES library for the rest of this work.

4.4 Kinematic results

In Figure 3, we present the high-resolution SINFONI kinematic maps of the central $3\times 3$ arcsecs of the galaxies resulting from the pPXF fits. The first column shows the S/rN map which we derived from the comparison between the pPXF fit and the input spectra (after applying the Voronoi binning). It visualizes well the quality of the pPXF fit and the quality of the data, as the S/rN is directly related to the errors of the kinematics, being large in the center and monotonically decreasing with radius. The S/rN maps show that our kinematics extraction works well (S/rN ¿ 30) within 1 arcsec which is the region that we used for our dynamical modeling. The next four columns show the velocity, velocity dispersion, $h_{3}$ and $h_{4}$ maps for each of our galaxies.

As expected from our selection criteria, the derived kinematics show mainly regular features. For each of our six galaxies, a clear rotation pattern is visible with maximal relative velocities ranging from 50 to 180 km/s (after subtracting the systemic velocity). The velocity dispersions show various patterns for the different galaxies. NGC 2784 and NGC 4281 contain a clear sigma increase within the isophotal center. The sigma peak in NGC 2784 has a size of about 0.3″  and goes up to 275 km/s, while we find a larger sigma peak in NGC 4281 ($\sigma\approx 310$ km/s). In NGC 3640 another clear feature is apparent: a slightly asymmetric velocity dispersion decrease in the center (down to 175 km/s) which spans the complete central region ($r<0.7\leavevmode\nobreak\ \arcsec$). This dip velocity is consistent with early work by Prugniel et al. (1988) and Davies et al. (1987). Prugniel et al. (1988) also point out that this galaxy might be in an advanced merger state which would significantly affect our dynamical models. Large scale signatures of this merger (such as shells) are also visible in the MATLAS images from Duc et al. (2015), also shown in Bonfini et al. (2018). However, Krajnović et al. (2011) analyzed the ATLAS3D kinematics of NGC 3640 with Kinemetry (Krajnović et al. 2006) and found only very small residuals and a very regular shape within one Re, indicating that the center of NGC 3640 is relaxed now. We, therefore, believe that our MBH mass measurement is robust and not likely affected by the advanced merger state (Prugniel et al. 1988). The velocity dispersion map of NGC 584 shows an hourglass-shape which can be attributed to a dynamically cold disc component (low velocity dispersions). Also, NGC 4570 shows signatures of a central disk. Its velocity dispersion goes up to 230 km/s, and we see maximal rotational velocities of 60 km/s which is fully consistent with the HST/Faint Object Spectrograph kinematics from van den Bosch et al. (1998). The velocity dispersion map of NGC 7049 is very unusual: it shows a very flat velocity dispersion profile without a clear sigma rise being visible in the kinematics of this galaxy.

The $h_{3}$ Gauss-Hermite moment maps show the typical anti-correlation to the velocity for each galaxy. The $h_{3}$ map of NGC 3640 may look chaotic at first glance, but the anti-correlation trend is also slightly visible here.

The visual comparison of the near-infrared central kinematic maps with the optical large-scale maps (Appendix: Fig. 15, 16) shows globally consistent results and similar trends even though we probe both very different scales and very different wavelength regions. The kinematic details of the SINFONI maps are generally not present on the large-scale kinematic maps, as the spatial resolution of the latter is comparable to the SINFONI FOV. In a second more quantitative comparison, we compared the Gauss-Hermite profiles from the four-moment pPXF fit of the two data sets. For the ”point-symmetric” velocity dispersion and $h_{4}$ moment, we averaged the bins within concentric circular annuli around the kinematic center and repeated this process with growing radius. The bins were chosen such that the luminosity-weighted center was within the respective annulus. We present the comparison in Figure 4. For some cases, we had to slightly shift the velocity dispersion and $h_{4}$ values of the SINFONI data (values are shown in Figure 4) in order to perfectly match the large and small-scale data. The shifts are about 5 % for three of our galaxies, no shifts for the remaining galaxies. Even before the shifts, all measured SINFONI velocity dispersions and most $h_{4}$ profiles were in very good agreement with the large-scale data. Some discrepancy can be seen in the $h_{4}$ profile of NGC 4281 which has a positive gradient for SINFONI and constant value for ATLAS${}^{\textrm{3D}}$. We believe that this discrepancy arises from the ATLAS${}^{\textrm{3D}}$ spatial resolution which flattens out the central features of the $h_{4}$ moment. Krajnović et al. (2018b) test the significance of the shifting with respect to the measured black hole mass and conclude that they add an uncertainty of about 80% of the measured black hole mass by shifting the velocity dispersion by about 8%. This means that we possibly add an uncertainty of 50% in mass for NGC 3640, NGC 4281 and NGC 7049.

5 DYNAMICAL MODELLING

We derived the central black hole masses of our target galaxies using two different and independent dynamical modeling methods: Jeans Anisotropic Models (JAM, Cappellari 2008) for constraining the parameter space and three-integral Schwarzschild (1979) orbital superposition models for deriving the final black hole masses. In the past, the Schwarzschild method has successfully been used to reproduce detailed models for spherical, axisymmetric and triaxial nearby galaxies. On the other hand, the JAM method is less general than orbit-based methods but far less computationally time-consuming. Furthermore, it provides a good description of galaxies based on two-dimensional stellar kinematics. Previous works, on almost forty galaxies, have shown that, although starting from different assumptions, both techniques provide generally consistent SMBH mass results (Cappellari et al. 2010; Seth et al. 2014; Drehmer et al. 2015; Thater et al. 2017; Feldmeier-Krause et al. 2016; Krajnović et al. 2018b; Ahn et al. 2018), such that modeling the observed stellar kinematics with both independent methods provides a more robust measurement. Recently, Leung et al. (2018) compared the results from both Schwarzschild and JAM models against circular velocities derived from molecular gas for 54 galaxies with CALIFA (Sánchez et al. 2012) integral-field stellar kinematics. They found that JAM and Schwarzschild recover consistent mass profiles (their Fig. D1). Moreover, JAM was found to recover more reliable circular velocities than the Schwarzschild models, especially at large radii where the gas velocities are more accurate (their fig. 8). Their study illustrates the fact that the reduced generality of the JAM method, with respect to the Schwarzschild method, is not necessarily a weakness and highlight the usefulness of comparing both methods as we do here.

5.1 The mass model

The gravitational potential of the galaxy is a composition of the potential of the stars, the potential from the central black hole which is assumed as a point mass, and the potential of dark matter. In order to find the mass of the central black hole, it is crucial to determine the stellar and dark matter contribution of the total galaxy mass as precisely as possible. The stellar mass density of the galaxy can be inferred from the galaxy luminosity density multiplied by the stellar M/L, which itself can be derived by modeling the stellar surface brightness of the galaxy. An efficient tool to provide an analytical description of the surface brightness of galaxies is the Multi-Gaussian Expansion (MGE) developed by Emsellem et al. (1994) and Cappellari (2002) in which a sum of two-dimensional concentric Gaussians parametrizes the galaxy surface brightness.

We performed the MGE modelling on both highly-resolved HST and deep wide-field ground-based SDSS (presented by Scott et al. (2013)) or CGS (Ho et al. 2011) imaging data, simultaneously, using the MgeFit Python package131313See footnote 2 Version 5.0 of Cappellari (2002). Except for NGC 7049, the SDSS images were in the r-band, while from HST we chose images taken with the WFPC2 camera in bands which matched the SDSS r-band best. NGC 7049 was only observed with the ACS camera in the F814W filter, and we matched it with I-band data from the CGS survey. We aligned the surface brightness profiles by re-scaling the large FOV imaging data to the central HST profiles and used the HST imaging for the photometric calibration. Foreground stars and nearby galaxies were carefully masked before applying this procedure. Furthermore, we had to apply a dust-correction to NGC 4281 and NGC 7049 to improve the modeling of the underlying galaxy surface brightness. Dust can have a significant effect on the stellar mass model as it alters the shape of the stellar surface brightness and dilutes the observed galaxy light. A careful dust correction is necessary to optimize the reproducibility from the model and the actual shape of the galaxy. We used the same method as in Cappellari et al. (2002) and Scott et al. (2013) to dust-correct the SDSS and CGS images and the dust-masking method outlined in Thater et al. (2017) to mask dust-rings visible in the HST images, which had only a single image available (for details see Appendix A). We also visually inspected the HST images of our galaxies for nuclear star clusters, but could not find any evidence. This is expected as galaxies with a mass of more than $10^{11}M_{\odot}$ usually do not harbor nuclear star clusters (Ferrarese et al. 2006; Wehner & Harris 2006; Seth et al. 2008; Graham & Spitler 2009).

The final MGE fits converge for between 10-12 Gaussian components centered on the galaxy nucleus and with the major axis aligned with the galaxy photometric major axis. For most of our lenticular galaxies, we can see a clear trend of the axial-ratio change with radius. They show rather round isophotes in the central bulge region and flattened and discy isophotes for larger radii due to the outer disc.

We converted the MGE parameters from pixel-space into physical units of L*⊙* pc*-2* following the guideline given by the MGE readme and Holtzman et al. (1995). For the transformation we needed to account for the absolute Vega magnitude of the sun141414http://mips.as.arizona.edu/~cnaw/sun.html M${}_{\rm F555W}=4.85$, M${}_{\rm F606W}=4.66$ and M${}_{\rm F814W}=4.15$. Furthermore, we corrected for the foreground Galactic Extinction applying the values found in the NASA/IPAC Extragalactic Database151515https://ned.ipac.caltech.edu/ which were derived by Schlafly & Finkbeiner (2011). The final MGE parameters are presented in Table 22: for each galaxy, we list the index of the Gaussian component, the surface brightness in units of L*⊙* pc*-2*, the Gaussian dispersion $\sigma_{j}$ in arcseconds and the axial ratios qj. In Figure 5, we show a comparison of our resulting best-fit MGE models and the observed HST WFPC2 and ACS images. Except for nuclear dust patterns (NGC 2784, NGC 4281, NGC 7049), the modeled MGE surface brightness are in good agreement with the surface brightness of each of the six galaxies. Especially for NGC 4281 a large dust mask had to be applied to correct the MGE model for the large amount of dust in this galaxy.

The next step for determining the mass model is the de-projection of the surface brightness into a three-dimensional luminosity density. Therefore, it is necessary to impose an assumption on the structure of the potential. For each of our target galaxies we adopted the assumption of an axisymmetric potential, such that, assuming a given inclination ($i>0$), the luminosity density can directly be deprojected from the MGE model. We used the built-in MGE regularization in order to bias the axial ratio of the flattest Gaussian to be as large as possible to prevent strong variations in the mass density of the MGE model. The MGE deprojection assumption does not remove the intrinsic degeneracy of the deprojection (Rybicki 1987; Gerhard & Binney 1996), which, especially at low inclination, can lead to major uncertainties and constitutes a fundamental limitation to the accuracy of any dynamical model (e.g. Lablanche et al. 2012). In the center of the galaxies, which is probed by our data, stars mainly contribute to the mean potential of the galaxy. This means that the galaxy mass density $\rho$ can simply be described as the product of the galaxy luminosity density and a dynamical mass-to-light ratio M/L. The gravitational potential generated by this mass density can then be obtained with the Poisson equation, $\nabla^{2}\Phi=4\pi G\rho$, and is one of the ingredients for the dynamical models in the next two sections. For further details regarding the MGE de-projection, we refer to the original work by Emsellem et al. (1994) and Cappellari (2002).

5.2 Jeans Anisotropic Models

The motion of a collection of stars in a gravitational field can be described by the Jeans (1922) equations. They provide the basis for the JAM method (Cappellari 2008), which predicts the second velocity moment by solving the Jeans and Poisson equations for the mass density derived from the MGE model. Projected along the line-of-sight of the model the second velocity moment is a function of four free parameters: the mass of the black hole $M_{BH}$, the anisotropy parameter $\beta_{z}$, the mass-to-light ratio M/L and the inclination angle i. The anisotropy parameter describes the orbital distribution by relating the velocity dispersion parallel to the rotation axis and in the radial direction: $\beta_{z}=1-\sigma_{z}^{2}/\sigma_{R}^{2}$ assuming that the velocity ellipsoid is aligned with cylindrical coordinates. We used the JAM method in order to model the second velocity moment in the potential defined by our MGE models, which is assumed to be axisymmetric. The modeled second velocity moment was then compared to the observed $V_{rms}=\sqrt{V^{2}+\sigma^{2}}$ with V being the mean velocity and $\sigma$ the velocity dispersion which was measured from the high-resolution SINFONI stellar kinematics (assuming a parametrization of the LOSVD of a simple Gaussian). Unlike the Schwarzschild models (Section 5.3), we here only fit the innermost high-resolution SINFONI kinematics to be robust against possible gradients in the M/L or the anisotropy.

We found the posterior distributions and the best-fitting values of the JAM parameters by applying a Bayesian framework in the form of Markov chain Monte Carlo (MCMC) inference method (Hastings 1970). We used the emcee software package Foreman-Mackey et al. (2013) which is a python implementation of the Goodman & Weare (2010) affine invariant Markov chain Monte Carlo ensemble sampler161616http://dfm.io/emcee/current/. JAM is generally fit to the data using Bayesian approaches and MCMC as this makes it easy to detect degeneracies between parameters and marginalize over uninteresting parameters (e.g., Cappellari et al. 2012; Barnabè et al. 2012; Cappellari et al. 2013b; Watkins et al. 2013; Cappellari et al. 2015; Mitzkus et al. 2016; Poci et al. 2016; Kalinova et al. 2017; Li et al. 2017; Bellstedt et al. 2018; Leung et al. 2018), and in the context of massive black hole determination by Krajnović et al. (2018b) and Ahn et al. (2018). For our dynamical JAM modeling, we followed a similar approach as Cappellari et al. (2013a). In the burn-in phase, a set of walkers explores the pre-defined parameter space, where each successive step is evaluated based on the likelihood of each walker. We used 100 walkers and tracked them for 200 steps until the fit converged. After the exploration of the parameter space, we continued the MCMC for 500 steps (post-burn-phase) and used the final walker positions to generate posterior distributions and model properties.

We built models with the four free parameters ($\log M_{BH}$, $\beta_{z}$,M/L,$i$) and compared them with the observed $V_{rms}$ using a $\chi^{2}$ statistic. The logarithmic likelihood probability of our data is defined as

[TABLE]

which is a sum over all good spaxels and where $\delta V_{rms}$ are the errors derived by the Monte Carlo simulations of the kinematic data and error propagation. In order to ensure that the fitting converges, we set reasonable priors on the parameters. We used uninformative priors (assumption of maximal ignorance) for the different parameters, which are uniform within the bounds of the likelihood function: $\log M_{BH}\in[4.8,9.8]$, $\beta_{z}\in[-1,+1]$, M/L $\in[0.1,20]$ and the inclination was allowed to vary over the full physical range (only limited by the flattening parameter $q_{\rm min}$ of the flattest Gaussian of the MGE model $\cos^{2}i=q^{2}$). We made sure that the MCMC chain converges by visually checking our burn-in plots and running several Markov chains.

In Appendix E (Figure 17), we present the MCMC posterior probability distributions of the different JAM model parameters for each galaxy. The contour plots show the projected two-dimensional distributions for each parameter combination and the histograms show the one-dimensional distributions for each parameter. As clearly indicated by the contour plots, our $M_{BH}$ and $\beta_{z}$ parameters are not degenerate for NGC 584, NGC 2784, NGC 3640 and NGC 4570, which shows that these measurements are robust. The $\beta_{z}$ parameters and the inclinations are naturally correlated but do not affect the black hole mass measurement. Generally, the derived inclinations are not well constrained and tend to be larger than expected from the literature. This is expected behavior as we only fit the central kinematics of our galaxies. That is why we decided to use the literature inclinations for the Schwarzschild modeling analysis. Furthermore, NGC 4281 and NGC 7049 show a degeneracy but still clearly constrain the black hole mass. We used the posterior probability distributions to calculate the best-fit value and their corresponding 3$\sigma$ uncertainties. The median values of the posterior distribution are given in Table 17. A visual comparison between the observed $V_{rms}$ profiles and the best-fit jam models of the MCMC routine are presented in Figure 6. In all cases the derived models reproduce the central peak of the observed $V_{rms}$ well, while the outer kinematics often suffer from scatter. Our derived best-fitting JAM MBH masses and M/L were used as an initial guess to constrain the Schwarzschild models.

5.3 Schwarzschild models

In our second dynamical modeling approach, we used the axisymmetric Schwarzschild code which was optimized for two-dimensional IFU data and described in Cappellari et al. (2006). The method is based on the numerical orbit-superposition method originally invented by Schwarzschild (1979) and further developed to fit stellar kinematics (Richstone & Tremaine 1988; Rix et al. 1997; van der Marel et al. 1998; Cretton et al. 1999). The basic idea of the Schwarzschild method is that the mass distribution of the galaxy is well described by the sum of time-averaged orbits in a stationary galaxy potential. The method consists of basically two steps which are repeated for each modeled black hole mass, respectively. First, assuming a stationary galaxy potential, a representative orbital-library is constructed from the galaxy potential which itself is derived from the mass density from Section 5.1. Regular orbits in axisymmetric potentials are characterized by three integrals of motion: the binding energy E, the vertical component of the angular momentum $L_{z}$ and a nonclassical third integral $I3$ introduced by Ollongren (1962), see also Richstone (1982) , which are equally sampled by the orbit library. We typically trace each orbit for 200 oscillations through the system to have a representative characteristic within the entire equilibrium phase of the galaxy.

In a second step, each orbit is projected into the plane of the observables and the complete set of orbits is combined to match the light distribution and the LOSVD of the galaxy by assigning a weight in a non-negative least-squared (NNLS) fit (Lawson & Hanson 1974). Compared to the JAM models, where we approximated the velocity second moments as the dispersion of a Gaussian, the Schwarzschild modeling method fits the full LOSVD.

We constructed the Schwarzschild models along a grid of radially constant dynamical mass-to-light ratio (M/L) and the mass of the black hole $M_{BH}$. We began the modeling procedure by running coarse parameter grids centered on the best-fitting parameters ($M_{BH}$, M/L) derived from the JAM models in Section 5.2. These models were improved iteratively by running finer and finer grids centered on the respective $\chi^{2}$ minimum of the coarse grid. Our final grids were then built with 21 $M_{BH}$ and 21 M/L equally spaced values for each galaxy. We only had to compute the orbit libraries for the different black hole masses as the orbits depend on the shape of the galaxy potential. The different M/L values only scale the potential and thus the orbit libraries can be re-scaled to match the different M/L a posteriori. Each orbit library consists of $21\times 8\times 7\times 2$ orbit bundles, which are composed of $6^{3}$ dithers, making in total 508 032 orbits per black hole mass. These orbit libraries were then fitted to the symmetrized stellar kinematics and to the photometric model in a NNLS fit and $\chi^{2}$ values were calculated by fitting our Schwarzschild models to both small and large-scale kinematics. We excluded the large-scale kinematics in the central $0.8\arcsec$ such that in the central regions only the more reliable high-resolution data was fitted. For the nnls fitting, we applied a regularization of $\Delta=10$ (analogous to Krajnović et al. (2009); van der Marel et al. (1998)) to impose an additional smoothing on the distribution function of the orbit weights. We present our final grids of Schwarzschild models for each of our six galaxies in Figure 7. Plotted on the grid is the $\chi^{2}$ distribution as a function of $M_{BH}$ and dynamical M/L from which we deduced the best-fitting parameters within $3\sigma$ significance ($\Delta\chi=11.8$). In order to smooth the topology of the $\chi^{2}$ contours, we applied the local regression smoothing algorithm LOESS Cleveland (1979), adapted for two dimensions (Cleveland & Devlin 1988) as implemented by Cappellari et al. (2013a, see footnote 3).

For each galaxy, we can constrain the upper and lower limit of the black hole masses. The best-fit values are presented in Table 17. Figure 7 also includes our JAM black hole mass measurements (MBH values within 99.7% intervals from posterior, namely 3$\sigma$) as grey shaded regions and the lowest possible black hole measurement based on the data resolution in combination with the sphere-of-influence argument (blue dashed line). NGC 3640, NGC 4281 and NGC 7049 have a clear overlap between the $3\sigma$ uncertainties of the JAM and Schwarzschild models meaning they are fully consistent with each other. For the remaining galaxies we measure slightly smaller black hole masses than with the JAM method. We note that the presented uncertainties on our black hole mass measurements are predominantly formal random errors from the dynamical modeling and as such they under-estimate the fuller systematic uncertainties which we discuss in Section 6.1. In Figure 8, we compare the $V_{RMS}$ maps between the SINFONI data and the Schwarzschild models for the best-fitting, a lower and higher MBH mass (just outside the $3\sigma$ contours) as well as the profiles along the x-axis. The different models are even visually very different, such that we can clearly constrain the upper and lower limit of the black hole mass. A full comparison between our observed (symmetrized) kinematic maps and the best-fitting Schwarzschild models with all our LOSVD parameters for both the SINFONI and large-scale data is shown in Appendix F. The models can reproduce all of the kinematic features very well, both on the high-resolution SINFONI data and the large-scale data.

NGC 4281 and NGC 7049 have an unusually large M/L, but roughly comparing the derived value of NGC 4281 (F606W-band) with the value from Cappellari et al. (2013a) who derived a value of 9.1 for the r-band by applying dynamical JAM models on the ATLAS${}^{\textrm{3D}}$ data only, our value is fully consistent.

6 DISCUSSION

The results that we recovered from our dynamical models are only robust when the assumptions on the models are valid. Therefore, we further investigated a number of systematic error sources that could have affected our results. In that respect, the choice of distance D does not influence our conclusions but sets the scale of our models in physical units. Specifically, lengths and masses are proportional to D, while M/L scales as D*-1*.

6.1 Systematic uncertainties

6.1.1 Variations in stellar populations

Various radial gradients have been found for different stellar population properties in early-type galaxies. For instance, early-type galaxies typically show color gradients, the central regions being redder than the galaxy outskirts (Peletier et al. 1990; Wu et al. 2005). Metallicities often follow a negative trend with radius (i.e., the metallicity decreases when the radius increases), while the age gradient is moderately flat (Kuntschner et al. 2010; Li et al. 2017). The mentioned gradients imprint their signature on the stellar M/L which is thus expected to increase towards the center. Furthermore, variations in the stellar IMF corresponding to a larger fraction of low-mass stars can have an additional effect on the M/L variation. Negative stellar M/L gradients were observationally confirmed for local early-type galaxies (e.g., recently in Boardman et al. 2017; Sarzi et al. 2018; Vaughan et al. 2018). Possibly problematic, in the previous section, we assumed the M/L to be constant for simplicity. However, ignoring the stellar M/L gradients can lead to overestimating the dynamical M/L and therefore also the central black hole mass (McConnell et al. 2013; Krajnović et al. 2018a). On the other hand, the stellar M/L usually runs contrary to the dark matter content which is low in the center but increases towards the outskirts of the galaxy. Therefore, including a nearly constant dynamical M/L must not always be a bad assumption in dynamical modeling (e.g. Thater et al. 2017) in particular, when modeling the stellar kinematics observed over a wide range of radial scales.

In our first dynamical modeling attempt, we assumed a constant dynamical M/L for both Jeans and Schwarzschild dynamical models. Comparing the dynamical M/L derived from the JAM models (where we only used the central kinematics ¡1.5″) and the Schwarzschild models, we noticed a significant (¿10%) difference for half of our sample: NGC 2784, NGC 4281 and NGC 4570. We considered that the dynamical M/L difference could be caused by stellar population variations. In order to study the effect of spatial variations in the stellar populations, we followed the same method as in McDermid et al. (2015) and Thater et al. (2017) and applied a mass-weighted stellar population synthesis for NGC 4570. We chose NGC 4570 for this test as its data had the best S/N and it did not suffer from dust contamination.

The ATLAS${}^{\textrm{3D}}$ IFU spectra of NGC 4570 were co-added in growth curves with increasing circular aperture sizes having radii between 0.5 and 25 arcsec and then fitted with a linear combination of MILES simple stellar population (SSP) model spectra (Vazdekis et al. 2010) using the pPXF routine. We used two different sets of template model spectra assuming either a unimodal initial mass function (IMF) of slope 1.30 (which equals a Salpeter (1955) IMF) or a Kroupa (2001) revised IMF. For each IMF choice, we used 350 SSP template spectra spanning a grid of 50 ages logarithmically spaced between 0.06 to 17.78 Gyr and 7 metallicities [Z/H] = [-2.32, -1.71, -1.31, -0.71, -0.40, 0.00, 0.22]. In addition, we also kept track of the stellar and stellar remnant mass $M_{*}$ and the r-band luminosity $L_{r}$ of each stellar model of the template library. Each of the template SSP spectra is assigned weights in the pPXF fit, which are smoothed out for models having similar ages and metallicities to ensure a smooth star formation history solution and suppress the noise in the final weights distribution. The smoothing is applied by adding a linear regularization to the pPXF fit which is chosen such that the difference in $\chi^{2}$ between regularised and non-regularised fit equals $\sqrt[]{2N}$, where N is the number of good pixels in the spectrum. We then calculated the mass-weighted stellar M/L for each radial bin using the tracked stellar mass and r-band luminosity from the SSP models and using equation (5) from Thater et al. (2017).

The derived M/L profiles of NGC 4570 for the two different IMFs, the metallicity and age profiles are shown in Figure 9. Within the effective radius, a negative M/L gradient in the order of 10-20 % of the central M/L is clearly visible, which has to be accounted for in the dynamical models. The gradient is very strong between 3″and 10″and flattens out for larger distances. Furthermore, while the shape of the M/L profile does not depend on the choice of the IMF, we note that their values differ by about $0.66\,M_{\odot}/L_{\odot}$ due to the ratio of high-mass to low-mass stars within the different IMFs. Recent papers suggest a trend in IMF with $\sigma_{e}$, in such a way that low-$\sigma_{e}$ (¡250 km/s) typically follow Kroupa-like IMFs, while galaxies with large $\sigma_{e}$ follow Salpeter-like or heavier IMFs (e.g. Cappellari et al. 2012, 2013a; Posacki et al. 2015; Li et al. 2017). Having velocity dispersions between 170 km/s and 245 km/s our sample galaxies thus would likely follow a Kroupa-like IMF but are located in the transition zone. In addition, IMF gradients have been found to follow the radial trend of the stellar metallicity (e.g. Martín-Navarro et al. 2015; Sarzi et al. 2018) which gives even more reasons to consider also bottom-heavy IMF forms. That is why we decided to derive the M/L for both IMFs and test them in our dynamical models.

The derived M/L values were then multiplied with the luminosity model MGE at the respective distance from the galaxy center (assuming an aperture size of the order of MGE $\sigma_{j}$ from Section 5.1) and included in the dynamical Schwarzschild models as mass density. We emphasize that we only included the M/L gradients in the Schwarzschild models as the JAM models only trace the galaxy potential within 1.5″where the stellar M/L is approximately constant. However, when constructing the Schwarzschild models, we also include stellar orbits from greater distances which could feel the effect of the observed M/L gradient. In order to account for the two possible IMFs, we ran the Schwarzschild grid for the two M/L profiles independently. We present the final black hole masses derived from the Schwarzschild models in combination with a variable M/L in Table 18. We find two main results from this analysis: 1) both IMFs give very consistent results which was expected as the shape of their M/L gradient is very similar, 2) including M/L variations in the Schwarzschild models reduces the derived black hole mass by about a factor of 1.5 (30%). The mass of the SMBH is decreased as more mass is included in the stellar component, and the impact on the black hole mass may have been even more important if we could have accounted for stellar population gradients down to the resolution of our SINFONI data. Our test agrees with McConnell et al. (2013) who have noticed that the MBH mass decreases by about 20-30% by taking M/L gradients into account. On the other hand, Cappellari et al. (2002) only found negligible variations when allowing for M/L gradients which were within the statistical uncertainties. While this test provides an interesting implication on the SMBH scaling relations, we will postpone a more detailed discussion for a future paper in the series when we can apply a detailed test to all 18 galaxies of the sample. This test will be crucial for the three galaxies of our sample that contain nuclear disks which are often accreted and thus likely have a different stellar populations and varying M/L gradients. Furthermore, it will be interesting to test if our dusty galaxies follow positive M/L gradients due to ongoing star formation and how much these gradients will affect the derived black hole measurements. A solution to the uncertainties introduced by possible unknown population gradients consists of allowing the total mass profile to differ from the distribution of the tracer population producing the kinematics (as done, e.g. in Mitzkus et al. 2016; Poci et al. 2016; Li et al. 2017). In black hole studies, this is generally done by allowing for a dark matter profile in addition to a stellar component (see Sec. 6.1.2), but the very same approach will account for gradients in the stellar M/L as well.

6.1.2 Dark matter

Our dynamical models only work under the assumption of self-consistence (mass follows light). Breaking this assumption by having significant amounts of dark matter in the center can lead to systematic changes in the black hole mass. (Gebhardt & Thomas 2009; Schulze & Gebhardt 2011; Rusli et al. 2013). We tested the significance of the dark matter in the central regions of our galaxies using the radial acceleration relation (McGaugh et al. 2016; Lelli et al. 2017). As long as the galaxies stay in the linear regime of the radial acceleration relation (gdyn ¿ g${}_{\mathrm{\mathrm{crit}}}=1.2\times 10^{-10}$ m s*-2*) it is expected that the dark matter does not contribute extensively to the galactic potential. The total acceleration can be derived from the gravitational potential by g${}_{\mathrm{dyn}}(R)=-\nabla\Phi_{\mathrm{tot}}(R)=V_{c}^{2}/R$ where $V_{c}$ is the circular velocity. We used the mass density (derived in Section 5.1) assuming the dynamical M/L of the best-fitting Schwarzschild model (Table 17) to calculate a model circular velocity at a radius equal to the edge of our large-scale kinematical data for each of our target galaxies. Our analysis yielded total accelerations between $8.8\times 10^{-9}$ and $9\times 10^{-10}$ m s*-2* with the smallest acceleration found in NGC 3640. All values lie well above the critical acceleration, and we conclude that our galaxies have likely a negligible contribution of dark matter in the central region which will not affect the dynamical modeling significantly. This is consistent with more direct estimates of the DM content of our galaxies from Cappellari et al. (2013a) and Poci et al. (2016). Furthermore, the total accelerations determined in our galaxies are consistent with the accelerations of other ATLAS3D early-type galaxies analyzed by Lelli et al. (2017).

6.2 Black hole - host galaxy scaling relations

We populated the M${}_{\mathrm{BH}}-\sigma_{\mathrm{e}}$ diagram with the compilation of dynamical black hole masses from Saglia et al. (2016) and Krajnović et al. (2018b). We then added our derived Schwarzschild MBH measurements in combination with the bulge effective velocity dispersions from Cappellari et al. (2013a)(see Table 1). The diagram is shown in Figure 10. Our measurements are located in the intermediate mass regime for early-type galaxies where the scatter is very tight. In Figure 10, we furthermore show the scaling relations derived in McConnell & Ma (2013), Saglia et al. (2016) and Savorgnan & Graham (2015). All of our black hole mass measurements follow the black hole scaling relation closely. Except for NGC 4281, our measurements are slightly below the scaling relation (but within the $1\sigma$ scatter of the relation), NGC 7049 deviating slightly more from the scaling relation.

The massive black holes in NGC 584 and NGC 3640 have already been measured indirectly in the literature. In their study Dullo & Graham (2014) recognized that NGC 584 and NGC 3640 show signatures for a partially depleted core which can be translated into a black hole mass of $M_{\mathrm{BH,dep}}=(1.95\pm 1.1)\times 10^{8}M_{\odot}$ for NGC 584 and $M_{\mathrm{BH,dep}}=(9.77\pm 1.1)\times 10^{7}M_{\odot}$ for NGC 3640. These values are consistent with our dynamical mass measurements of $M_{\mathrm{BH,dyn}}=(1.34\pm 0.49)\times 10^{8}M_{\odot}$ for NGC 584 and of $M_{\mathrm{BH,dyn}}=(7.73\pm 0.51)\times 10^{M}_{\odot}$ for NGC 3640. We also note that our dynamical mass measurements derived from the Jeans modeling match perfectly with the MBH masses derived from the depleted cores.

While our measurements follow the general trend of previous mass measurements, a systematic offset seems to emerge between MBHs in early-type and late-type host galaxies, the latter being significantly lower. Graham & Scott (2013) note that this offset is also seen for barred versus non-barred galaxies (barred galaxies having larger velocity dispersion). The authors note that most of the late-type galaxies on the M${}_{\mathrm{BH}}-\sigma_{\mathrm{e}}$ relation are actually barred, and it is not clear at this moment if the departure of the late-type galaxies from the scaling relation for the early-type galaxies is driven by bars or is typical for all late-type galaxies. Furthermore, black hole masses measured via $\rm H_{2}O$ megamasers (e.g., Greene et al. 2010) in possibly barred galaxies also seem to be systematically lower than dynamical black hole mass measurements (which is nicely visualized in Figure 1 of van den Bosch 2016, see also Davis et al. 2018a, b). In Section 6.1.1 we have recognized that by taking into account a variable M/L the dynamical mass measurements could shift down by a factor of about 1.5. The radial M/L variation might even be more important for late-type galaxies, except for cases where the MBH is estimated by directly observing tracers within the MBH sphere of influence, such as H2O megamasers. We will investigate this implication in a future paper of this series. Independent from the M/L variation, together with the recently published dynamical mass measurements by Krajnović et al. (2018b) our measurements strengthen the idea of early-type galaxies having more massive black holes than (barred) late-type galaxies and thus following different M${}_{\mathrm{BH}}-\sigma_{\mathrm{e}}$ relations (e.g., Terrazas et al. 2017; Davis et al. 2018a; Martín-Navarro et al. 2018).

7 SUMMARY AND CONCLUSION

In this work, we have presented the black hole mass measurement of six nearby early-type galaxies (NGC 584, NGC 2784, NGC 3640, NGC 4281, NGC 4570 and NGC 7049). Our measurements are based on AO-assisted K-band SINFONI IFU observations complemented by ground-based IFU data from MUSE, VIMOS, and SAURON from the ATLAS${}^{\textrm{3D}}$ survey. All of our target galaxies show regular rotation and except for NGC 3640 and NGC 7049 a velocity dispersion increase towards their center. NGC 3640 contains a velocity dispersion dip of about 30 km/s, while NGC 7049 seems to have a flat velocity dispersion profile. This finding is consistent with the kinematic features of the large-scale SAURON data from the ATLAS${}^{\textrm{3D}}$ survey and our VIMOS and MUSE data. We combined our kinematic results with photometric mass-models based on the composition of HST, SDSS and CGS survey data to build dynamical models to measure the mass of the central black holes. We constrained the parameter space of possible masses and mass-to-light ratios using Jeans Anisotropic Modelling on our SINFONI data and then created axisymmetric orbit-superposition Schwarzschild modeling based on both central and large-scale IFU data to derive robust results. We derive black hole masses of $(1.3\pm 0.5)\times 10^{8}M_{\sun}$ for NGC 584, $(1.0\pm 0.6)\times 10^{8}M_{\sun}$ for NGC 2784, $(7.7\pm 5)\times 10^{7}M_{\sun}$ for NGC 3640, $(5.4\pm 0.8)\times 10^{8}M_{\sun}$ for NGC 4281, $(6.8\pm 2.0)\times 10^{7}M_{\sun}$ for NGC 4570 and $(3.2\pm 0.8)\times 10^{8}M_{\sun}$ for NGC 7049 which fit well with the recent black hole - $\sigma_{\mathrm{e}}$ scaling relations.

For three galaxies we find a slight discrepancy in the derived dynamical M/L of the two different methods, the central values being larger than the M/L derived from the combination of small-scale and large-scale data. Dynamical models typically assume a constant M/L for simplicity reasons which are usually not the case in observed galaxies. To test this assumption, we derive the stellar M/L profile from stellar population modeling for the test case of NGC 4570, which does not suffer from dust contamination and has the best quality data. The stellar population modeling shows a negative gradient of about 20% within the effective radius of the galaxy which is based on variations in stellar age and metallicity. This negative gradient is then included in the dynamical Schwarzschild models, and we derive a black hole mass which is lower by almost 30%, irrespective of further possible stellar M/L re-scaling due to radially constant IMF variations. We conclude that the inclusion of M/L variations has an effect of the order of the general uncertainty of the measurement, but it should be included in dynamical models to lower the systematic uncertainties which are still very large in dynamical modeling. As was already suspected in different works, this has an interesting implication on the black hole scaling relations as different dynamical methods still suffer from partially inconsistent results. We caution that careful additional analysis of the effect of M/L variations on dynamical models is urgently needed in future studies.

Acknowledgements

The authors want to thank Aaron Barth of the Department of Physics and Astronomy, University of California, Irvine, for sharing the large scale imaging data of NGC 584, NGC 2784 and NGC 7049 from the CGS survey with us. MC acknowledges support from a Royal Society University Research Fellowship.

Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 075.B-0495(A), 078.B-0464(B), 079.B-0402(B), 080.B-0015(A), 085.B-0221(A), 091.B-0129(A), 097.A-0366(B), 291.B-5019(A). Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the Data Archive at 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 This research has made use of the NASA/IPAC Extragalactic Database which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.

Appendix A Dust correction and masking

As the presence of dust can have a crucial effect on the galaxy mass modeling, we had to correct and mask the dust polluted image pixels before constructing the mass models. NGC 4281 and NGC 7049 contain extended nuclear dust rings which are well visible in the HST images (see Figure 11). Furthermore, we found a small dust ring in the HST image of NGC 2784. As carefully tested in Thater et al. (2017) we applied a dust mask for the HST small scale images for all three galaxies and a dust correction for the SDSS large scale image of NGC 4281.

A.1 HST images

For NGC 2784, NGC 4281 and NGC7049 the nuclear dust ring reaches into the very central regions of the galaxies. It was therefore necessary to also correct for dust in the HST small scale images. The careful dust correction in the HST images is crucial as we probe the direct vicinity of the black hole with them. However, due to the lack of HST images in different bands, we applied the dust masking method which was developed in Thater et al. (2017). The method is based on the assumptions 1) that the presence of dust attenuates the light emitted within the galaxy and 2) the major part of the galaxy surface brightness is not significantly affected by dust. In order to derive the dust affected images we fitted the lower envelope of the characteristic surface brightness profile with an appropriate function (four parameter logistic function). Masked were all pixels which had a surface brightness below this fit. This method corrects patchy dust structures which can be clearly distinguished from the unattenuated regions, but misses thin dust structures and dust layers (see Figure 11).

A.2 Large scale images

In order to correct the dust attenuation in the SDSS and CGS images we applied a method that was developed by Carollo et al. (1997) and applied and further advanced by Cappellari et al. (2002), Krajnović et al. (2005) and Scott et al. (2013). The main assumption is that dust between the observer and the stellar emission can be folded into a screen which dims the observed light wavelength-dependent. Due to the extinction, the observer sees a change in colour in the dust-affected regions of the galaxy which are assumed to have the same intrinsic colour as their adjacent areas. Using the Galactic extinction law (Schlegel et al. 1998) we derived the r-band extinction for NGC 4281 from the colour excess between the g- and i-band images $A_{r}=1.15\,E(g-i)$ . In order to derive the colour excess we created the respective colour profile as function of the logarithm of the semi-major axis distance along ellipses with fixed position angle and ellipticity. The colour profile of NGC 4281 is presented in Figure 12 showing a very slight colour gradient with radius and the signatures of the nuclear dust ring in the reddened pixels above. We fitted a robust linear function to the colour gradient to reduce the influence of the dust-affected pixels. The colour excess $E\,(g/r-i)$ of each pixel is now computed by subtracting its assumed intrinsic colour value (approximated by the linear fit) from its measured value. All pixels above an arbitrary threshold value, here chosen to be $E\,(g-i)>0.02$ mag for NGC 4281, are assumed to have a significant dust extinction and were corrected using the Galactic extinction law. Figure 13 also shows quantitatively how much of the measured flux was corrected (where 0.1 means 10% (blue) and 0.25 means 25% (orange)). The largest correction was approximately 35% of the measured flux. Figure 12 also shows the colour profile of NGC 7049. Based on this plot and an additional visual check of the image we realized that the dust content in this galaxy is mostly concentrated in the centre (within the HST PC FoV) and we decided to only apply the dust-masking of the HST images.

Appendix B Determination of the HST, NIFS and VIMOS PSF

B.1 HST spatial resolution

We generated a PSF image for the HST WFPC2/PC and ACS images using the online available Tiny Tim HST Modelling tool191919http://www.stsci.edu/hst/observatory/focus/TinyTim (Krist & Hook 2001) taking into account the imaging filter, the central position of the galaxy on the CCD chip and assuming the spectrum of a K giant star. While we could use the HST WFPC2 images as they were provided by Tiny Tim, the procedure was slightly more complicated for the ACS images. The ACS camera has the disadvantage to be mounted away from the optical axis of HST which distorts the observed image significantly giving it the shape of a sheared rectangle (TinyTim manual202020http://tinytim.stsci.edu/static/tinytim.pdf). Distortion is not corrected by the internal optics, but must be accounted for during the data-reduction, both for the science image and the PSF image. In order to do so, we used a method applied by Rusli et al. (2013). We substituted the distorted PSF image into the flat-fielded galaxy ACS image, replacing the centre of the galaxy. We then corrected the full substituted ACS image using the drizzelpac task in Astroconda with the same input parameters as during the science image data reduction. We then cut out the PSF image from the resulting distortion-corrected ACS mosaic and applied the same following procedures as for the WFPC2 images.

The resulting PSF image was parametrized by a sum of normalized concentric circular Gaussians using the MGE method Cappellari (2002). All MGE parameters and relative weights of each Gaussian and for each galaxy are given in Table 7.

B.2 SINFONI and VIMOS spatial resolution

The spatial resolution of the IFU data sets a limit on the scales that can be probed with the central dynamics of our target galaxies. As such the resolution provides a quality argument of our black hole mass measurements which has to be carefully evaluated. When no point source is present in the FoV, the reconstructed IFU images need to be compared with reference images of significantly higher resolution. We used HST imaging data (used for the MGE models) in order to determine the PSF of the SINFONI observations (e.g. Shapiro et al. (2006); Davies (2008); Krajnović et al. (2009). Deconvolved HST images (a product from the MGE modelling) are convolved with a PSF (and such degraded) until they match the collapsed IFU image. A simple parametrisation of the SINFONI PSF can be obtained using a circular and concentric double Gaussian composed of a narrow and broad component of different relative weights. In order to compare the different images, they were aligned by rotating the images such that the major and minor axes would coincide with the vertical and horizontal image axes. Furthermore, after the convolution the MGE model was re-binned to the same pixel scale as the respective IFU observation. The best fitting parameters were found by minimising the residual between the convolved MGE model and the reconstructed SINFONI image using the python routine mpfit.

We furthermore performed a number of different tests with the PSF determination routines which involved changing the flux scales between the HST and collapsed IFU image, varying the rotation angle, the size of the fitted image and the kernel size. This led to additional systematic uncertainties of about 0.05 arcsec for the narrow Gaussian component. A comparison of the main axes between the SINFONI light profiles and convolved MGE models is presented in Fig. 14. Even though we probe very different wavelength regions, the fits agree remarkably well. The best-fitting parameters of the PSF fits are also given in Fig. 14 and summarized in Table 6.

B.3 Strehl ratio

The Strehl ratio is an indication of the effect of wavefront aberrations on the optical quality of the observations. It can be determined by calculating the ratio between the peak intensity of the measured PSF and the peak intensity of the ideal diffraction limited PSF assuming an ideal working LGS AO. We obtained the FWHM of the narrow Gaussian component of the SINFONI observation from the PSF determination (see Table 6). The diffraction limit of the telescope can be calculated as $\theta=1.22\times\lambda/D$ where $\lambda$ specifies the wavelength of the observed light and $D$ is the diameter of the primary mirror. At 2.3 microns, the diffraction limit of the VLT telescope (D= 8.2 m) is approximately 0.07 ″. We created normalized 2D Gaussians with the resolution of the SINFONI observation and the diffraction limit. Dividing the peak intensities resulted into Strehl ratios of around 10 % for our SINFONI observations (see Table 6).

Appendix C MGE parametrisation of target galaxies

Appendix D Large scale kinematics

Appendix E JAM models

Appendix F Comparison of the Schwarzschild dynamical models with the symmetrised data