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

(1)

Astronomy& Astrophysics manuscript no. axi_paper ESO 2019c October 5, 2019

Six new supermassive black hole mass determinations from

adaptive-optics assisted SINFONI observations

Sabine Thater

1

, Davor Krajnovi´c

1

, Michele Cappellari

2

, Timothy A. Davis

3

, P. Tim de Zeeuw

4, 5

,

Richard M. McDermid

6

_{, and Marc Sarzi}

7, 8

1 _{Leibniz-Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany,e-mail: sthater@aip.de}

2 _{Sub-Department of Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK}

3 _{School of Physics & Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK}

4 _{Sterrewacht Leiden, Leiden University, Postbus 9513, 2300 CA Leiden, The Netherlands}

5 _{Max Planck Institute for Extraterrestrial Physics (MPE), Giessenbachstrasse 1, D-85748 Garching b. München, Germany}

6 _{Department of Physics and Astronomy, Macquarie University, Sydney, NSW 2109, Australia}

7 _{Centre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane,}

Hatfield, Hertfordshire AL10 9AB, UK

8 _{Armagh Observatory and Planetarium, College Hill, Armagh, BT61 9DG, UK}

Received ... 2018/ Accepted ?

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 (σ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 of (1.3 ± 0.5) × 108_M

for NGC 584, (1.0 ± 0.6) × 108M

for NGC 2784, (7.7 ± 5) × 107_M

for NGC 3640, (5.4 ± 0.8) × 108Mfor NGC 4281, (6.8 ± 2.0) × 107Mfor NGC 4570 and

(3.2 ± 0.8) × 108_M

for NGC 7049 at 3σ confidence level, which are consistent with recent MBH- σescaling 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.

Key words. galaxies: individual: NGC 584, NGC 2784, NGC 3640, NGC 4281, NGC 4570, NGC 7049 – galaxies: kinematics and dynamics – galaxies: supermassive black holes

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 dynami-cal 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 σ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-σ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´c 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´c et al. (2018a) and Mezcua et al. (2018) recently reported an up-bending of the scaling relations

(2)

Table 1. The sample

Galaxy Type Distance Linear scale MK Re σe σ0 log(Mbulge) i Large Scale (Mpc) (pc arcsec−1) (mag) (arcsec) (km s−1) (km s−1) log(M) (◦)

(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) NGC 584 S0 19.1 ± 1.0 93 −24.19 16.5 191 ± 5 216 10.48 51 MUSE NGC 2784 S0 9.6 ± 1.8 47 −23.31 20.1 201 ± 6 243 10.44 66 VIMOS NGC 3640 E3 26.3 ± 1.7 128 −24.60 38.5 176 ± 8 173 11.00 68 SAURON NGC 4281 S0 24.4 ± 2.2 118 −24.01 24.5 227 ± 11 314 10.88 71 SAURON NGC 4570 S0 17.1 ± 1.3 83 −23.48 17.9 170 ± 8 209 10.18 88 SAURON NGC 7049 S0 29.9 ± 2.6 145 −25.00 17.7 245 ± 8 266 11.02 42 VIMOS

Notes – Column 1: Galaxy name. Column 2: Morphological type (de Vaucouleurs et al. 1991). NGC 584 had been misclassified in the earlier work and we adopt the classification by Huang et al. (2013) here. Column 3: Distance to the galaxy (taken from Cappellari et al. (2011) for SAURON/ATLAS3Dgalaxies or NED for VIMOS and MUSE galaxies), the uncertainties were calculated by dividing the NED standard deviations by √Nwhere N is the number of measurements. Column 4: Linear scale derived from the distance. Column 5: 2MASS total K-band magnitude (Jarrett et al. 2000). Column 6: Effective radius derived from B-band (CGS) or r-band (ATLAS3D_{) imaging data. The values were taken from Ho et al. (2009) or Cappellari et al. (2013a) for ATLAS}3D_galaxies. Column 7: Effective velocity dispersion derived by co-adding the spectra of the large-scale optical IFU data in elliptical annuli of the size of the effective radius. For ATLAS3D_{galaxies taken from Cappellari et al. (2013a). Column 8: Central velocity dispersion} derived by co-adding the spectra of the high-spatial-resolution SINFONI IFU data in elliptical annuli within one arcsec. Column 9: Bulge mass calculated by multiplying the bulge-to-total ratios from Krajnovi´c et al. (2013) for ATLAS3D_{galaxies or Gao et al.} (2018) for the remaining galaxies with the total dynamical mass from Cappellari et al. (2013a) or this paper. Column 10: Inclination from Cappellari et al. (2013a) for ATLAS3D_{galaxies and Ho et al. (2011) for remaining galaxies. For NGC 584 from Laurikainen} et al. (2010). Column 11: Large scale kinematics data which is used for the Schwarzschild dynamical models. The SAURON data comes from the ATLAS3Dgalaxy survey.

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 MBHby 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 (≈ 300 km/s) end. First results, based on observations with NIFS and GEM-INI, were published in Thater et al. (2016) and Krajnović et al. (2018b). Unlike many other MBHmeasurements 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 MBHmeasurements, especially if the sphere of influence is hardly resolved.

(3)

2. The Sample

The six galaxies analyzed in this paper belong to our SMASH-ING galaxy sample to dynamically determine black hole masses in the nearby universe. Three of our target galaxies were selected from the ATLAS3D _{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 re-maining galaxies were observed with the VIMOS or MUSE in-struments. 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 MBHrange. The main properties of our six sample galaxies are summarized in Table 1.

Our target galaxies were selected based on a number of require-ments for a successful MBHdetermination. An important crite-rion 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 de-pends on the mass of the black hole MBHand the velocity disper-sion of the galaxy within an effective radius σeand is defined as rSoI= GMBH/σ2ewhere G is the gravitational constant. We cal-culated an estimated value for rSoIusing black hole masses based on the MBH−σerelation from Tremaine et al. (2002)1 and the ATLAS3D _{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´c et al. 2009; Thater et al. 2016). With the goal to gain the best possible resolution, we utilized the AO mode from the SINFONI instru-ment 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 galax-ies 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 un-suitable 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 infor-mation of the central galaxy region to constrain the wide range of different stellar orbit families and large-scale IFU data to con-strain the global galaxy characteristics are essential for a pre-cise measurement. The IFU data is complemented by imaging data from HST and ground-based telescopes to construct a de-tailed mass model of the host galaxy. In the following section, we present the different observations from the IFUs towards the imaging data.

1 _{The 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

se-lect galaxies that were most likely to provide robust MBH 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.

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 Ob-servations 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 In-frared Faint Field Imaging (SPIFFI) assisted by the adaptive op-tics (AO) module, Multi-Application Curvature Adaptive Opop-tics (MACAO) (Eisenhauer et al. 2003; Bonnet et al. 2004). We ob-served each galaxy at K-band grating (1.94 - 2.45 µm) providing a spectral resolution of R ∼ 4000 and a pixel scaling of 100 mas leading to a total field-of-view (FOV) on the sky of about 3.2 × 3.200 _{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 high-est spatial resolution possible. In the ideal case, the LGS mode still requires a natural guide star to correct for the tip-tilt dis-turbances 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 pro-vided a slight improvement to the natural seeing. Our observa-tions show typical Strehl ratios of about 10 % (see Table 3). The observations were performed using the object-sky-object nod-ding 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 re-construct the data cubes of the individual observations. This sci-ence 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 de-scribed in the SINFONI instrument handbook. The observa-tions were reduced using the ESO SINFONI pipeline (ver-sion 2.4.8, Modigliani et al. 2007) in combination with addi-tional 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

(4)

Table 2. Details of the SINFONI observing runs

Galaxy Date PID Pixel scale N N Texp 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

Notes. Column 1: Galaxy name. Column 2: Dates of the observations. Column 3: Identification number of the Proposal. Column 4: Spatial pixel scaling of the observation. Column 5: Number of available single exposure frames. Column 6: Number of single exposure frames in the combined data cube. Column 7: Combined exposure time in hours. Column 8: Adaptive optics mode applied for the data, either using a natural guide star (NGS) or a laser guide star (LGS).

science cubes.

For the telluric correction of the near-infrared spectra, we wrote a Python script to apply the same method as described in Kra-jnovi´c 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 tar-get. The telluric stars were either solar-like G0-2V stars or hotter B2-5V stars in an unsystematic order. We used the Python ver-sion 6.06 of the penalized Pixel fitting software2_(pPXF, Cappel-lari & Emsellem 2004) as upgraded in CappelCappel-lari (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 so-lar template (Livingston & Wallace 1991)3 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 λcorrected = (1+ vhelio/c) × λuncorrectedwhere c is the speed of light and vhelio is the projected heliocentric velocity which was calculated from the ESO pipeline for each data frame. The heliocentric correc-tion 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 ob-serving 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

2 _{http://purl.org/cappellari/software}

3 _{ftp://nsokp.nso.edu/pub/atlas/photatl/}

galaxy would not always fall onto the same pixels of the detec-tor and thus adding systematic uncertainties. It was possible to identify the center of the galaxy with an accuracy of about one pixel (=0.0500) 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 see-ing 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´c et al. (2009) and Thater et al. (2016). The algo-rithm defines a new square pixel grid and interpolates the science frame to this grid. Flux values of the final data cubes were cal-culated as median flux values of the single data frames. Finally, we obtained 3 × 300data cubes with 0.0500pixel scale for the 100 mas SINFONI observations.

(5)

Fig. 1. Example of the spectral resolution inhomogeneity across the SINFONI detector. The spectral resolution for each spaxel was derived from arc line observations of NGC 584. The spectral resolution varies significantly in the vertical direction of the detector with values ranging from 5.5 Å to 7.7 Å FWHM.

of 6.8 Å FWHM (λ/∆λ = 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 method4 _{(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 es-tablish typical bin sizes of< 0.100 in the center, while 0.3-0.400 diameter for bins at a radius larger than 100_.

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 3.

4 _{See footnote 2}

Table 3. SINFONI spatial resolution

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

Notes. The SINFONI PSF of the data was parametrized by a double Gaussian with a narrow and broad component. The parameters are given in the following columns. Column 1: Galaxy name. Column 2: FWHM of the narrow Gaussian component. Column 3: FWHM of the broad Gaussian component. Column 4: Relative intensity of the narrow com-ponent. Column 5: Strehl ratio of the data.

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 900s on-source integrations each yielding a field of view of 60 × 6000(≈ two effective radii of NGC 584). In addi-tion, 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 ev-ery 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 func-tion for each slice. Each on-target observafunc-tion was reduced us-ing 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 correc-tion. The sky spectrum was applied to all three on-target frames, where we let the pipeline model the sky lines based on the in-put 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 pro-cedure. In the final cube each pixel has the size of 0.200×0.200 and a spectral sampling of 1.25 Å per pixel. For our MBH de-termination we only needed the high S/N central 3000_{× 30}00 _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.500in the center and 300at radii larger than 1200_.

3.2.2. VIMOS IFU data

(6)

The VIMOS data reduction was performed by Lagerholm et al. (2012) making use of the ESO pipeline5 _{(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 re-moval. As described in Lagerholm et al. (2012), they also cor-rected 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 VI-MOS data to a target S/N of 60, obtaining bin sizes of 0.500in the galaxy center and 2-300_{at radii larger than 7.5}00_.

3.2.3. SAURON IFU data

NGC 3640, NGC 4281 and NGC 4570 are part of the ATLAS3D 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 Willia4.2-m 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 3300× 4100with a sampling of 0.9400× 0.9400square pixels, cov-ering about 1-2 effective radii of our target galaxies. A detailed description of the stellar kinematics extraction of the ATLAS3D sample is given in Cappellari et al. (2011). The resulting velocity maps of NGC 3640, NGC 4281 and NGC 4570 were already pre-sented in Krajnovi´c et al. (2011), and we show the full kinematic set of these galaxies in Figure D.1. 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 ATLAS3Ddata6which 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 Ad-vanced Camera for Survey (ACS, Ford et al. 1998) data from the ESA Hubble Science Archive, which generates automatically re-duced 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 im-age was reprocessed by applying the drizzlePac7_{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 Sur-vey 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 ATLAS3D collaboration (Scott et al. 2013).

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

Notes. Column 1: Galaxy name. Column 2: programme identification number. Column 3 and 4: Camera on HST and the filters in which the data were taken.

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 method8 _{(pPXF, Cappellari & Emsellem 2004;} Cappellari 2017). pPXF fits the galaxy spectrum by con-volving 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 σ 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, σ) 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 statis-tical 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

5 _{http://www.eso.org/sci/software/pipelines/}

6 _{http://purl.org/atlas3d}

7 _{http://www.stsci.edu/hst/HST_overview/drizzlepac}

(7)

2.15 2.20 2.25 2.30 2.35 2.40 2.45 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 Normalized Flux

NGC 0584, S/rN=77

SINFONI

2.15 2.20 2.25 2.30 2.35 2.40 2.45 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

NGC 2784, S/rN=56

2.15 2.20 2.25 2.30 2.35 2.40 2.45 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

NGC 3640, S/rN=79

2.15 2.20 2.25 2.30 2.35 2.40 2.45 Wavelength [ m] 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 Normalized Flux

NGC 4281, S/rN=94

2.15 2.20 2.25 2.30 2.35 2.40 2.45 Wavelength [ m] 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

NGC 4570, S/rN=75

2.15 2.20 2.25 2.30 2.35 2.40 2.45 Wavelength [ m] 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

NGC 7049, S/rN=86

5000 5500 6000 6500 Wavelength [Å] 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 Normalized Flux

NGC 0584, S/rN=276

MUSE

4500 4750 5000 5250 5500 5750 6000 Wavelength [Å] 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

NGC 2784, S/rN=179

VIMOS

4500 4750 5000 5250 5500 5750 6000 Wavelength [Å] 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2

NGC 7049, S/rN=189

VIMOS

Fig. 2. Integrated SINFONI, MUSE and VIMOS spectra and pPXF fits of our target galaxies. The integrated spectra (black solid lines) were obtained by summing up all spectra of the IFU data cubes and fitted using the pPXF routine (red lines) in order to derive an optimal template. The fitting residual between spectrum and best fitting model are shown as green dots and are shifted up by 0.5 (0.6 for the bottom panels). Regions which were masked in the fit (often due to emission lines or insufficient sky subtraction) are indicated as grey shaded regions and their residuals are indicated in blue.

using χ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: ≈ 2.5 − 5 km s−1versus ≈ 5 − 10 km s−1_{and velocity dispersion: 2.5 − 6 km s}−1_versus 6 − 12 km s−1).

1.0 0.5 0.0 0.5 1.0

arcsec

0.15

0.00

0.15 NGC 4281

25

55

85 180 0

180 160 235 310

0.15

0.00

0.15

1.0 0.5 0.0 0.5 1.0

arcsec

0.15

0.00

0.15 NGC 4570

30

45

60 70 0

70 160 195 230

0.15

0.00

0.15

1.5 1.0 0.5 0.0 0.5 1.0 1.5

arcsec

0.15

0.00

0.15 NGC 7049

30

55

80 100 0

100 220 250 280

0.15

0.00

0.15

1.0 0.5 0.0 0.5 1.0

arcsec

0.15

0.00

0.15

Fig. 3. SINFONI stellar kinematics (derived from CO bandhead spectroscopy) of our target galaxies (from top to bottom) NGC 584, NGC 2784,

NGC 3640, NGC 4281, NGC 4570 and NGC 7049. From left to right the panels show maps of signal-to-residual noise (S/N), mean velocity (V),

velocity dispersion (σ) and the Gauss-Hermite moments h3and h4. The black contours indicate the galaxy surface brightness from the collapsed

(9)

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 µ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. Further-more, 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 ∼ 2.25 µ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 (σtemp ≈ 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 ∼ 2.2 µm as none of our stellar templates could match the line strength fully. As ? 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 ATLAS3Dkinematic extraction. We re-extracted the kinematics for VIMOS data as the Lager-holm 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 em-pirical spectra (MILES, Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011) stellar library9 _{(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 disper-sion of the MILES template library to be σ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

9 _{http://miles.iac.es/} 0 1 2 3 200 250

[k

m

/s]

NGC0584

0 1 2 3 0.1 0.0 0.1

h4

MUSE Sinfoni 0 1 2 3 200 250

m

/s]

NGC4570

0 1 2 3 0.1 0.0 0.1

h4

ATLAS3d Sinfoni -0.05 0 1 2 3

Radial bins [arcsec]

250 275

[k

m

/s]

×1.05

NGC7049

0 1 2 3

Radial bins [arcsec]

0.1 0.0 0.1

h4

VIMOS Sinfoni 0.06

Fig. 4. Comparison of velocity dispersion and h4 profiles for the

SIN-FONI (red) and the respective large scale data (blue). The values were averaged within circular annuli around the kinematic centre. The error range of the averaged values in the radial bins are calculated via error propagation and are shown as shaded regions. Applied shifts in the SIN-FONI maps are marked by the values in the upper right corner of each panel.

(10)

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, how-ever, discern a difference in the extracted values between the two stellar templates with the MILES velocity dispersions be-ing systematically lower (10-20 km/s) and thus, more consis-tent 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. Compar-ing 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 × 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, h3 and h4 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.300 _{and goes up to} 275 km/s, while we find a larger sigma peak in NGC 4281 (σ ≈ 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.700). 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´c et al. (2011) analyzed the ATLAS3D kinematics of NGC 3640 with Kinemetry (Krajnovi´c 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 h3 Gauss-Hermite moment maps show the typical anti-correlation to the velocity for each galaxy. The h3map 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 kine-matic maps with the optical large-scale maps (Appendix: Fig. D.1, D.2) 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 com-parable 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 h4 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 h4values 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 h4profiles were in very good agreement with the large-scale data. Some discrepancy can be seen in the h4profile of NGC 4281 which has a positive gradient for SINFONI and constant value for ATLAS3D. We believe that this discrepancy arises from the ATLAS3D _spatial resolution which flattens out the central features of the h4 moment. Krajnovi´c 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

(11)

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 accu-rate (their fig. 8). Their study illustaccu-rates the fact that the reduced generality of the JAM method, with respect to the Schwarzschild method, is not necessarily a weakness and highlight the useful-ness 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 package10 _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. (2016) 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 1011Musually 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.

10 _{See footnote 2}

We converted the MGE parameters from pixel-space into physical units of Lpc−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 sun11_M

F555W = 4.85, MF606W = 4.66 and MF814W = 4.15. Fur-thermore, we corrected for the foreground Galactic Extinction applying the values found in the NASA/IPAC Extragalactic Database12 which were derived by Schlafly & Finkbeiner (2011). The final MGE parameters are presented in Table C.1: for each galaxy, we list the index of the Gaussian component, the surface brightness in units of Lpc−2, the Gaussian dispersion σjin 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 ρ 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, ∇2Φ = 4πGρ, 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 MBH, the anisotropy parameter β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: βz= 1−σ2z/σ2Rassuming that the velocity ellipsoid is aligned with cylindrical coordinates. We used the

11 _{http://mips.as.arizona.edu/~cnaw/sun.html}

NGC 7049 (ACS)

Fig. 5. Isophotal maps of the WFPC2 and ACS images of our target galaxies within a FoV of 20 × 2000

. In the bottom right of each panel we show

a cutout of the central 3 × 300

. The contours of our best-fitting MGE model (red) are superimposed on the HST images. For the models, foreground stars and close galaxies were masked. For NGC 2784, NGC 4281 and NGC 7049 a dust-correction of the internal dust rings had to be applied

before MGE modelling their surface brightness. The MGE models were build from the combined photometric information of HST (r< 1000_{) and}

wide-field of view data (r> 1000

) from ATLAS3D_{and the CGS (Ho et al. 2011) survey.}

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 Vrms =

√

V2+ σ2 _{with V being} the mean velocity and σ the velocity dispersion which was measured from the high-resolution SINFONI stellar kinematics (assuming a parametrization of the LOSVD of a simple Gaus-sian). 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 frame-work 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 invari-ant Markov chain Monte Carlo ensemble sampler13. 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

13 _{http://dfm.io/emcee/current/}

et al. 2018; Leung et al. 2018), and in the context of massive black hole determination by Krajnovi´c 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.

(13)

1.5 1.0 0.5 0.0 0.5 1.0 1.5 160 180 200 220 240 260 280 300 NGC 584 Data JAM 1.5 1.0 0.5 0.0 0.5 1.0 1.5 160 180 200 220 240 260 280 300 Vrms (k m /s) 1.5 1.0 0.5 0.0 0.5 1.0 1.5 220 240 260 280 300 320 NGC 2784 1.5 1.0 0.5 0.0 0.5 1.0 1.5 220 240 260 280 300 320 Vrms (k m /s) 1.5 1.0 0.5 0.0 0.5 1.0 1.5 140 160 180 200 220 240 NGC 3640 1.5 1.0 0.5 0.0 0.5 1.0 1.5 140 160 180 200 220 240 Vrms (k m /s) 1.5 1.0 0.5 0.0 0.5 1.0 1.5 200 250 300 350 400 NGC 4281 1.5 1.0 0.5 0.0 0.5 1.0 1.5 200 250 300 350 400 Vrms (k m /s) 1.5 1.0 0.5 0.0 0.5 1.0 1.5 170 180 190 200 210 220 230 240 250 260 NGC 4570 1.5 1.0 0.5 0.0 0.5 1.0 1.5 170 180 190 200 210 220 230 240 250 260 Vrms (k m /s) 1.5 1.0 0.5 0.0 0.5 1.0 1.5 arcsec 220 230 240 250 260 270 280 290 300 NGC 7049 1.5 1.0 0.5 0.0 0.5 1.0 1.5 arcsec 220 230 240 250 260 270 280 290 300 Vrms (k m /s)

Fig. 6. Comparison of the Vrms profiles between the SINFONI data

(blue) and best-fitting JAM models (green) along the major (left) and minor (right) axis. The green shaded region shows JAM models with varying black hole mass by a factor of 1.3 either larger or smaller than the best-fitting mass.

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 MBH∈ [4.8, 9.8], βz∈ [−1,+1], M/L ∈ [0.1, 20] and the inclination was allowed to vary over the full physical range (only limited by the flattening parameter qminof the flattest Gaussian of the MGE model cos2_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 E.1), we present the MCMC pos-terior probability distributions of the different JAM model pa-rameters 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 MBHand βz parameters are not degenerate for NGC 584, NGC 2784, NGC 3640 and NGC 4570, which shows that these mea-surements are robust. The βzparameters and the inclinations are naturally correlated but do not affect the black hole mass mea-surement. Generally, the derived inclinations are not well con-strained 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 incli-nations 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 correspond-ing 3σ uncertainties. The median values of the posterior distri-bution are given in Table 5. A visual comparison between the observed Vrmsprofiles and the best-fit jam models of the MCMC routine are presented in Figure 6. In all cases the derived mod-els reproduce the central peak of the observed Vrmswell, 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 Lz 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.

(14)

0.5 1.0 1.5 2.0 2.5 3.0

1e8

4.5

5.0

5.5

6.0 M/

L (

M

¯

/L

¯

)

11.8 23.6 47.2 94.4 188.8 188.8 377.6 377.6 755.2 755.2

NGC 584

0.5 1.0 1.5 2.0

1e8

5.5

6.0

6.5

7.0

7.5

8.0

11.8 23.6 47.2 94.4 94.4 188.8 188.8 377.6

NGC 2784

0.0 0.5 1.0 1.5 2.0 2.5

1e8

3.6

3.8

4.0

4.2

4.4

4.6

4.8

11.8 11.8 _23.6 47.2 47.2 94.4 94.4 188.8 188.8

NGC 3640

0.2 0.4 0.6 0.8 1.0

M

BH

(M

¯

)

1e9

M

BH

(M

¯

)

1e8

10.5

11.0

11.5

12.0

12.5

13.0

13.5

11.8 23.6 23.6 47.2 47.2 94.4 94.4 188.8 377.6

NGC 7049

Fig. 7. Grids of Schwarzschild models (indicated by the black dots) over different mass-to-light ratios M/L and black hole masses MBH. The

best-fitting model, derived as the minimum of χ2_{, is indicated by a large red circle. Contours are the}_∆χ2_{= χ}2₋_χ2

minlevels where the thick green

contour shows the 3σ level of the two-dimensional distribution. In addition, we have added the 3-sigma limits on the best-fitting black hole masses of the JAM models (grey shaded regions). The dashed blue line indicates the mass of the black hole which has the radius of the sphere of influence of half the resolution of our LGS AO data (inferred from the narrow component of the AO PSF), which is approximately the lowest black hole measurement that we expect to be detectable based on our data.

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 MBH. We began the modeling procedure by running coarse parameter grids centered on the best-fitting parameters (MBH, 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 χ2 _{minimum of} the coarse grid. Our final grids were then built with 21 MBHand 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

(15)

1

0

1 NGC 0584

M¯

200

220

240

260 V

RM S

[k

m

/s]

1

0

1 NGC 2784

1 .

03

×

10

8

M¯

3 .

19

×

10

7

M¯

1 .

61

×

10

8

M¯

220

240

260

280 V

RM S

[k

m

/s]

1

0

1 NGC 3640

7 .

73

×

10

7

M¯

0 M¯

1 .

88

×

10

8

M¯

160

180

200 V

RM S

[k

m

/s]

1

0

M¯

250

300

350 V

RM S

[k

m

/s]

1

0

M¯

180

200

220 V

RM S

[k

m

/s]

1

0

1

0

/s]

Fig. 8. Comparison of the VRMS=

√

V2_{+ σ}2_{maps from the SINFONI data and the Schwarzschild models. Each row shows the maps of one galaxy,}

respectively. From left to right we present the observed symmetrized VRMSfrom the SINFONI data, and the VRMSmaps of the Schwarzschild models

from the best-fitting, a too low and a too high MBHas well as the profiles along the x=0 axis. The too low (blue) and too high (orange) black hole

masses are chosen to be just outside of the 3σχ2_{contours. All models are shown at the respective best-fitting M/L. The high and low mass models}

are clearly ruled out for all galaxies.

(1979), adapted for two dimensions (Cleveland & Devlin 1988) as implemented by Cappellari et al. (2013a, see footnote 3).

(16)

Table 5. Summary of dynamical modelling results

JAM Schwarzschild

Galaxy MBH M/L β i χ2/DOF MBH M/L χ2/DOF rSoI/σPSF

(×108_M ) (M/L) (◦) (×108M) (M/L) (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) NGC 584 1.93±0.06 _5.4±0.1 _0.05±0.01 ₈₉±3 _4.35 _1.34±0.49 _5.4±0.2 _0.99 _1.6 NGC 2784 1.69±0.1 _7.7±0.3 _0.04±0.04 ₇₇±13 _1.86 _1.03±0.54 _6.7±0.7 _1.08 _1.8 NGC 3640 0.99±0.13 _4.2±0.1 _0.14±0.04 ₈₅±9 _13.90 _0.77±0.51 _4.2±0.2 _1.65 _1.1 NGC 4281 4.91±0.15 _12.5±0.2 _−0.03±0.02 ₇₅±9 _9.93 _5.42±0.80 _9.3±0.3 _3.98 _2.1 NGC 4570 1.19±0.09 6.6±0.1 0.22±0.02 74±1 3.98 0.68±0.20 5.5±0.1 1.87 1.1 NGC 7049 3.16±0.84 11.4±0.4 0.04±0.04 44±10 3.18 3.17±0.84 11.9±0.3 1.38 1.6

Notes. Column 1: Galaxy name. Column 2-6: Parameters of the JAM models (black hole mass, mass-to-light ratio, velocity anisotropy parameter and stellar mass of the galaxy). Column 7-9: Parameters of the Schwarzschild models (black hole mass and mass-to-light ratio in the HST band

specified in Table 4). Column 10: Comparison of the black hole sphere-of-influence (calculated with the central velocity dispersion σ0) and the

spatial resolution of the observations (measured by the narrow component of the AO PSF).

Table 5. Figure 7 also includes our JAM black hole mass mea-surements (MBH values within 99.7% intervals from posterior, namely 3σ) as grey shaded regions and the lowest possible black hole measurement based on the data resolution in combina-tion with the sphere-of-influence argument (blue dashed line). NGC 3640, NGC 4281 and NGC 7049 have a clear overlap be-tween the 3σ uncertainties of the JAM and Schwarzschild mod-els meaning they are fully consistent with each other. For the re-maining galaxies we measure slightly smaller black hole masses than with the JAM method. We note that the presented uncer-tainties 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 VRMS maps between the SINFONI data and the Schwarzschild models for the best-fitting, a lower and higher MBH mass (just outside the 3σ contours) as well as the profiles along the x-axis. The di ffer-ent 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 ATLAS3D_{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´c 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. 2016) in particular, when modeling the stellar kinematics observed over a wide range of radial scales.