WISDOM Project - I: Black Hole Mass Measurement Using Molecular Gas Kinematics in NGC 3665

TL;DR

This study measures the supermassive black hole mass in galaxy NGC 3665 using molecular gas kinematics from millimeter-wave observations, providing precise estimates and confirming consistency with galaxy property correlations.

Contribution

First to model SMBH mass in NGC 3665 using 3D molecular gas kinematics with high-resolution millimeter observations.

Findings

SMBH mass estimated at approximately 5.75 x 10^8 solar masses.

Molecular gas kinematics consistent with previous methods and galaxy property correlations.

Systematic uncertainties quantified for SMBH mass and stellar mass-to-light ratio.

Abstract

As a part of the mm-Wave Interferometric Survey of Dark Object Masses (WISDOM) project, we present an estimate of the mass of the supermassive black hole (SMBH) in the nearby fast-rotator early-type galaxy NGC 3665. We obtained Combined Array for Research in Millimeter Astronomy (CARMA) B and C array observations of the $^{12}$CO$(J=2-1)$ emission line with a combined angular resolution of $0".59$. We analysed and modelled the three-dimensional molecular gas kinematics, obtaining a best-fit SMBH mass $M_{\rm BH}=5.75^{+1.49}_{-1.18} \times 10^{8}$ $M_{\odot}$, a mass-to-light ratio at $H$-band $(M/L)_{H}=1.45\pm0.04$ $(M/L)_{\odot, H}$, and other parameters describing the geometry of the molecular gas disc (statistical errors, all at $3\sigma$ confidence). We estimate the systematic uncertainties on the stellar $M/L$ to be $\approx0.2$ $(M/L)_{\odot, H}$, and on the SMBH mass to be $\approx0.4\times10^{8}$ $M_{\odot}$. The measured SMBH mass is consistent with that estimated from the latest correlations with galaxy properties. Following our older works, we also analysed and modelled the kinematics using only the major-axis position-velocity diagram, and conclude that the two methods are consistent.

Tables

Table 1. Table 1: NGC 3665 properties.

Parameter	Value	Reference
Morphology	fast rotator early type	1
Position		2
RA (J2000.0)	${11}^h{24}^m{43}^s.624$
DEC (J2000.0)	$38°45\prime 46\stackrel{\prime \prime }{.}278$
Systemic velocity (km s-1)	2069	3
Position angle ($°$)	26	4
Inclination angle ($°$)	${69.9}_{-0.39}^{+0.51}$	4
Distance (Mpc)	$34.7\pm 6.8$	5,6
Linear scale (pc arcsec-1)	$167\pm 33$	5

Table 2. Table 2: CARMA observation parameters.

Parameter	B array	C array
Date	Nov-Dec 2013	Apr 2013
On-source time (min)	$1610.6$	$290$
Phase center:
RA(J2000.0)	${11}^h{24}^m{43}^s.6$
DEC(J2000.0)	$38°45\prime 46\stackrel{\prime \prime }{.}278$
Primary beam	$31″$
	LSB	USB
Frequency coverage (GHz)	$214.404$–$215.373$	$228.576$–$229.545$
Velocity resolution (km s-1)	$10.9$	$10.2$

Table 3. Table 3: MGE components of the HST NICMOS F160W and 2MASS H 𝐻 H -band images.

$j$	$I_j$2222	${\sigma }_j$	$q_j$
	($L_{\odot ,H}$ pc-2)	(arcsec)
1……….	$25551.5$22	$0.227$	$0.515$
2……….	$21118.8$22	$0.661$	$0.608$
3……….	$7436.97$	$1.31$	$0.887$
4……….	$12016.7$22	$2.17$	$0.576$
5……….	$5862.67$	$4.76$	$0.837$
6……….	$741.344$	$11.3$22	$0.441$
7……….	$807.669$	$19.2$22	$0.780$
8……….	$212.118$	$48.6$22	$0.821$

Table 4. Table 4: Model parameters.

Parameter	Search Range	Best Fit	Error ($1\sigma$ conf.)	Error ($3\sigma$ conf.)
SMBH mass (${10}^8{M}_{\odot }$)	$0.01$–$50.12$	$5.75$	$+0.42,-0.38$	$+1.49,-1.18$
Stellar $M/L$ ($M/L_{\odot ,H}$)	$0.10$–$4.00$	$1.45$	$\pm 0.01$	$\pm 0.04$
Molecular gas disc:
Centre X offset (arcsec)	$-3.50$–$3.50$	$-0.25$	$\pm 0.01$	$\pm 0.02$
Centre Y offset (arcsec)	$-3.50$–$3.50$	$0.02$	$\pm 0.01$	$\pm 0.03$
Inclination ($°$)	$67.00$–$89.00$	$69.90$	$\pm 0.20$	$\pm 0.61$
Position angle ($°$)	$0$–$50$	$26$	$\pm 0$	$\pm 1$
Centre velocity offset (km s-1)	$-50.00$–$50.00$	$-8.12$	$\pm 0.75$	$\pm 2.50$
Velocity dispersion (km s-1)	$1.00$–$20.00$	$12.53$	$\pm 0.74$	$\pm 2.09$
Luminosity scaling	$10.00$–$200.00$	$86.16$	$\pm 2.91$	$\pm 8.63$
Scale length (arcsec)	$1.00$–$7.00$	$2.11$	$\pm 0.01$	$\pm 0.03$
Central void radius (arcsec)	$0.01$–$0.90$	$0.38$	$\pm 0.03$	$\pm 0.04$

Full text

††volume: in press††pubyear: 2017††pagerange: WISDOM Project – I: Black Hole Mass Measurement Using Molecular Gas Kinematics in NGC 3665–WISDOM Project – I: Black Hole Mass Measurement Using Molecular Gas Kinematics in NGC 3665

1 Introduction

Substantial improvements in observational capabilities over the past two decades have allowed to measure the mass of supermassive black holes (SMBHs) in nearby galaxies using a few different methods across a range of wavelengths, revealing that SMBHs are ubiquitous and most likely present in every stellar spheroid (see, e.g., Kormendy & Ho 2013 for a review). Black holes are now recognized to play a fundamental but as yet not fully understood role in the growth and evolution of galaxies.

As the angular resolution of observations increases, so does the number of reliable SMBH mass measurements, and several correlations have now emerged between SMBH mass and host galaxy properties such as bulge mass (e.g. Magorrian et al., 1998; Marconi & Hunt, 2003; Beifiori et al., 2012; Kormendy & Ho, 2013), total luminosity (e.g. Läsker et al., 2014) and bulge velocity dispersion (e.g. Ferrarese & Merritt, 2000; Gültekin et al., 2009; McConnell & Ma, 2013; Kormendy & Ho, 2013). This last so-called $M_{\rm BH}-\sigma$ relation is the tightest, but all empirical correlations suggest that galaxy evolution and SMBH growth are closely connected. Numerical simulations reproduce these correlations and provide suggestive evidence that this co-evolution likely involves self-regulation mechanisms (e.g. AGN feedback; Silk & Rees 1998; Di Matteo et al. 2008; Barai et al. 2014).

Another fact supporting the co-evolution paradigm is that as the number of dynamically-measured SMBH masses increases, different relationships for different types of galaxies are starting to emerge. For example, McConnell & Ma (2013) found that different relations best fit samples of early- and late-type galaxies, the SMBH mass being typically two times larger at a given velocity dispersion for early-type galaxies. The current empirical correlations are however based on only $\approx 80$ objects, not enough to reliably divide objects and contrast the relations of different galaxy samples. This shortcoming is compounded by the fact that almost $70\%$ of these objects are early-type galaxies, a problem arising from the very limited number of methods available to measure SMBH masses, and the fact that each method is biased toward objects with certain characteristics. We discuss this problem below, and the need for more measuring techniques applicable across a broader range of galaxy types.

Dynamical SMBH mass measurements have so far relied on only a handful of methods. Stellar dynamical measurements have the obvious advantage that every galaxy has a substantial population of stars moving exclusively under the influence of gravity, and stars are easily probed at optical and near-infrared (NIR) wavelengths. The modelling methods are well established and applicable to many galaxies; $\gtrsim 60\%$ of current SMBH mass measurements rely on stellar dynamics (see, e.g., Dressler & Richstone, 1988; van der Marel & van den Bosch, 1998; Cappellari et al., 2002; Gebhardt et al., 2003; Valluri et al., 2005; McConnell et al., 2011; van den Bosch et al., 2012; Rusli et al., 2013). Nevertheless, the dynamical modelling tools are relatively complex and often restricted to axisymmetric objects, dust can easily perturb measurements, and relatively high spectral resolution is required to probe the higher order moments of the line-of-sight velocity distributions. These constraints result in samples heavily biased toward early-type galaxies.

Ionised gas dynamics is also used to measure SMBH masses, and is typically probed through nebular emission lines at optical/NIR wavelengths, that are present in most galaxies (see, e.g., Ferrarese et al., 1996; Macchetto et al., 1997; van der Marel & van den Bosch, 1998; Cappellari et al., 2002; Neumayer et al., 2007; de Francesco et al., 2008; Walsh et al., 2013). Challenges for this method however include non-gravitational forces (e.g. shocks) and significant turbulent motions superimposed on (quasi-)circular motion. Observations with multiple slits are often unable to fully characterise the kinematics and overcome these challenges, but significant improvements have been made with integral-field spectroscopic (IFS) observations (see, e.g., Neumayer et al., 2007; Walsh et al., 2016).

The dynamics of accretion discs containing megamasers offers another way forward to measure SMBH masses (see, e.g., Miyoshi et al., 1995; Kuo et al., 2011). The precision is nearly on par with that in our own Milky Way (e.g. Ghez et al., 2008), occasionally allowing to positively rule in SMBHs (as opposed to other potential compact dark objects), but the rarity of appropriate megamaser systems ($\sim 1\%$ of objects searched; see Lo 2005; van den Bosch et al. 2016) makes this at best a niche method. Another limitation of this method is that megamasers have only been detected in Seyfert 2 and low-ionisation nuclear emission-line region (LINER) galaxies so far, containing SMBHs of relatively low masses. The reason for this bias is unknown, but it again suggests a difficulty to build up large samples.

More recently, the dynamics of molecular gas probed at mm/sub-mm wavelengths has emerged as a very promising method (see, e.g., Davis et al., 2013; Onishi et al., 2015; Barth et al., 2016a, b), in particular because of the exquisite angular resolution and sensitivity afforded by the Atamaca Large Millimeter/sub-millimeter Array (ALMA). As for ionised gas, molecular gas can in principle be affected by non-gravitational forces and turbulence (e.g. García-Burillo et al., 2003; Smajić et al., 2015), although in practice (giant) molecular clouds tend to move ballistically over much of their orbits, and molecular gas observations are unaffected by dust. Of course, some galaxies have no molecular gas detected (as is the case for many early-type galaxies; Young et al. 2011), and others none in their very centre (Helfer et al., 2003; Alatalo et al., 2013), but this method has enormous potential. For instance, Davis (2014) showed that the SMBH mass of $\approx 35,000$ local galaxies could be measured with ALMA when it reaches full capability, and the sphere of influence (SOI; $R_{\rm SOI}\equiv GM_{\rm BH}/\sigma^{2}$, where $G$ is the gravitational constant) of the largest SMBHs ($M_{\rm BH}\geq 10^{8.5}$ $M_{\odot}$) is spatially resolvable across the whole of cosmic time. Needless to say, the molecular gas method also has the potential to redress the current bias against late-type galaxies in $M_{\rm BH}-\sigma$ studies (now only $\approx 30\%$ of the sample).

Both Davis et al. (2013) and Onishi et al. (2015) derived the SMBH mass without directly detecting the Keplerian motion at the galaxy centre. This was due to the particular molecular gas disc morphology (in the case of Davis et al. 2013) and severe beam smearing (in the case of Onishi et al. 2015). The SMBH masses were constrained to good accuracy (e.g. 20% error for Onishi et al. 2015), although higher-angular resolution would have allowed to further constrain the measured values. Barth et al. (2016a) detected the Keplerian upturn at the centre of the nearby galaxy NGC 1332, but the SMBH mass was not very well constrained due to the SOI resolved along the disc’s major axis but unresolved along its minor axis. Higher-angular resolution observation in Barth et al. (2016b) constrained the SMBH mass to $10\%$ accuracy.

Building on small pilot projects (Davis et al., 2013; Onishi et al., 2015), we have recently started the mm-Wave Interferometric Survey of Dark Object Masses (WISDOM) project. This project aims to benchmark and test the molecular gas dynamics method, develop tools and best practice, and exploit the growing power of ALMA to better populate and thus constrain SMBH – galaxy scaling relations. This paper is the first of a series, and introduces the tools and fitting procedures developed so far, in the context of estimating the SMBH mass in the nearby fast-rotator early-type galaxy NGC 3665. We use the method initially employed by Davis et al. (2013), but further extend it to exploit the full three-dimensional data cubes.

This paper is stuctured as follows. The galaxy selection, observations, and data reduction and analysis are described in Section 2. The SMBH mass measurement method is explained in Section 3. Section 4 contains a discussion of galaxy morphology, AGN activity, possible sources of error on the SMBH mass, and a comparison of the SMBH masses derived using different fitting methods. Our main results are summarized in Section 5.

2 Data

2.1 NGC 3665

NGC 3665 is a nearby fast-rotator early-type galaxy (Emsellem et al., 2011). We adopt a distance of $34.7$ Mpc, estimated from the Tully-Fisher relation by Theureau et al. (2007), yielding a scale of $\approx 167$ pc arcsec*-1*. This distance is in good agreement with the estimate of $33.1$ Mpc from its recession velocity111The choice of the distance does not influence our conclusions but sets the scale of our models in physical units. Specifically, lengths and masses scale as $D$, while $M/L$s scale as $D^{-1}$. by Cappellari et al. (2011). Basic properties of the galaxy are summarised in Table 1. Twin jets emanating from the galaxy nucleus were observed with arcsecond resolution at $1.5$ GHz using the Very Large Array (VLA; Parma et al. 1986), while a point-like structure was observed with $2$ milli-arcsecond resolution at $5$ GHz using the Very Long Baseline Array (VLBA; Liuzzo et al. 2009), both proving the existence of a central SMBH. We adopt the position of this point-like structure as the galaxy centre: RA=11:24:43.624, Dec=38:45:46.278 (Liuzzo et al., 2009).

Hubble Space Telescope (HST) observations in $H$-band (NICMOS F160W; see the left panel of Figure 1) show prominent and regular dust lanes circling the galaxy centre, suggesting a gaseous disc-like structure extending to a radius of at least $8\arcsec$ from the nucleus. Previous CARMA observations revealed a regularly rotating molecular gas disc of the same radius (Alatalo et al., 2013). We therefore expect a central SMBH surrounded by a relaxed molecular gas disc, ideal to measure the SMBH mass using molecular gas kinematics.

NGC 3665 does not have a SMBH mass measurement yet, but has an effective stellar velocity dispersion of $\sigma_{\rm e}=216\pm 10$ km s*-1*, reliably measured through integral-field spectroscopy (Cappellari et al., 2013). Adopting the $M_{\rm BH}-\sigma$ relation of McConnell & Ma (2013), this suggest a SMBH mass of $\approx 3\times 10^{8}$ $M_{\odot}$ and thus a SMBH SOI of $\approx 30$ pc or $\approx 0\aas@@fstack{\prime\prime}2$, that is 3 times smaller than the synthesised beam size of the current observations (see Section 2.2).

2.2 Observations and Data Reduction

NGC 3665 was observed at the wavelength of the 12CO$(J=2-1)$ line using CARMA in the B (baselines $63$–$947$ m) and C (baselines $26$–$370$ m) arrays. Observations were carried out from April 11th to 21st 2013 (B array) and from November 26th to December 14th 2013 (C array). Total on-source time was 1610.6 min in the B array and 290 min in the C array. The receivers were tuned to cover the frequency range $214.404$–$215.373$ GHz and $228.576$–$229.545$ GHz in the lower (LSB) and upper (USB) sidebands, respectively, with $4$ spectral windows per range. With $31$ channel per spectral window, the frequency resolution was thus $7.812$ MHz per channel. The field of view (full width at half maximum, FWHM, of the primary beam) at these frequencies was $31\arcsec$ for the $10$-m antennae. The observational parameters are summarized in Table 2.

We followed the data reduction and analysis method described in Alatalo et al. (2013) using the Multichannel Image Reconstruction Image Analysis and Display (MIRIAD) package (Sault et al., 1995). The visibility data were first edited and calibrated using Mars, MWC349 and 3C273 as flux calibrators and 1153+495 as phase calibrator. The bandpass calibrator was 3C279 or 3C273. We then subtracted the continuum by assuming the CO emission to be present from $-450$ to $450$ km s*-1* with respect to the galaxy systemic velocity of $2069$ km s*-1* (Cappellari et al., 2011), and subtracting a linear fit to the line-free channels. The resulting cube was cleaned using the MIRIAD task MOSSDI2 with a threshold of 1.5 times the rms noise, measured in line-free regions of the cube. For imaging, we set the robustness parameter to $0.5$, yielding a synthesised beam FWHM of $0\aas@@fstack{\prime\prime}60\times 0\aas@@fstack{\prime\prime}56$ ($\approx 100$ pc $\times$ $93$ pc) at a position angle of $-81\degr$, that was properly sampled with $0\aas@@fstack{\prime\prime}2\times 0\aas@@fstack{\prime\prime}2$ pixels. This provides a compromise between angular resolution and sensitivity, and the rotational motion of the molecular gas in the galaxy nucleus is then clearly detected. We did not bin the velocity channels, yielding a velocity resolution of $7.8$ MHz or $\approx 10$ km s*-1* per channel, and the average rms noise per channel was $4.1$ mJy beam*-1*. Molecular gas emission was finally detected from $-370$ to $380$ km s*-1*.

Integrated intensity (moment 0) and intensity-weighted (mean) velocity (moment 1) maps can be created directly from this cube, but as most of the cube is devoid of emission, the resulting maps are of poor quality. Instead, we optimised the moments by first Hanning-smoothing the data cube in velocity and then smoothing it spatially with a Gaussian of FWHM equal to that of the synthesised beam. A mask was then created by selecting all pixels in the smoothed data cube above a threshold of $0.75$ times the rms noise in each channel. The adopted integrated intensity and mean velocity maps, shown in Figure 1, were then created by calculating the moments of the original unsmoothed data cube within the mask region only.

The integrated intensity map reveals a centrally-concentrated molecular gas distribution, rapidly decreasing with radius. However, the inner $\approx 2\arcsec$ clearly show two separate concentrations on either side of the nucleus (see Figure 1, left panel). This suggests a void in the very centre of the galaxy, where the intensity may remain above zero simply due to the angular extent of the synthesised beam. This central hole is confirmed by our modelling and further discussed in Section 4.1. The mean velocity map (Figure 1, right panel) reveals very regular disc-like rotation with a total velocity width of $\approx 750$ km s*-1* and no evidence of any significant non-circular motion, warp or kinematic twist.

Interestingly, the molecular gas distribution and dust lanes revealed by HST are not very well associated with each other (see Figure 1, left panel). The kinematic major axis of the molecular gas and the major axis of the dust lanes are however aligned, suggesting that the molecular gas disc and dust are nevertheless in the same plane. Both also align well with the large-scale photometric major axis, suggesting that this plane is also the large-scale equatorial plane of the galaxy. The differing distributions of the molecular gas and dust may thus be due to missing flux in the interferometric data, due to the lack of truly short baselines. Indeed, comparing the flux from our observations integrated over the CO$(2-1)$ Institut de Radioastronomie Millimetrique (IRAM) 30-m telescope beam ($42.3$ Jy km s*-1*) to an actual IRAM 30-m integrated flux measurement ($67.1$ Jy km s*-1*; Young et al. 2011), our CARMA high-resolution observations may be resolving out $\approx 40\%$ of the flux in extended structures.

3 Supermassive Black Hole Mass Estimate

In this section, we describe the procedures employed to measure the SMBH mass in NGC 3665 and state our results. To summarise, we modelled the three-dimensional stellar mass distribution of the galaxy by deprojecting a two-dimensional model of the observed surface brightness and assuming a constant mass-to-light ratio $M/L$. The circular velocity curve arising from this mass model and a putative SMBH was then fed into a code simulating the resulting data cube, taking into account the molecular gas distribution and instrumental effects. The SMBH mass was then determined by simply comparing a range of models to the observations.

3.1 Velocity Model

The galaxy mass distribution is assumed to be the sum of a central SMBH and the large-scale stellar body of the galaxy. The SMBH is treated as a point mass whose mass is free. For the stars, we combine HST (NICMOS F160W) and Two Micron All-Sky Survey (2MASS) $H$-band images, allowing to accurately trace the stellar surface brightness to a radius of $\approx 40\arcsec$. We adopt the Multi Gaussian Expansion (MGE) method (Emsellem et al., 1994; Cappellari et al., 2002) and fit this two-dimensional image with a sum of Gaussians (MGE_FIT_SECTORS procedure222Available from http://purl.org/cappellari/software of Cappellari (2002)). Given an inclination each Gaussian can be deprojected analytically, and the three-dimensional light distribution of the model can thus be trivially reconstructed. Here all the Gaussians are constrained to have the same position angle and inclination, resulting in an axisymmetric light model.

The point spread function (PSF) of each image is also fit with a sum of (circular) Gaussians and used as input during the MGE fit to obtain a deconvolved light model of the galaxy. We use the Tiny Tim package (version 6.3) developed by Krist et al. (2011) to measure the HST NICMOS PSF of the F160W filter. The FWHM of the 2MASS $H$-band PSF was assumed to be $2\aas@@fstack{\prime\prime}8$333http://spider.ipac.caltech.edu/staff/roc/2mass/seeing/seesum.html. The prominent dust lane seen in the HST image is masked to mitigate the effects of dust obscuration, and the region outside a radius of $6^{\prime\prime}$ in the HST image is ignored and overwritten with the 2MASS image. Figure 2 shows a comparison of the best-fitting MGE model and the observed surface brightness distribution of the galaxy, from both 2MASS and HST. We follow the photometric calculation described in Section 5.2 of Thatte et al. (2009) to convert the flux units from counts pixel*-1* second*-1* to $L_{\odot}$ pc*-2*, and the resulting MGE parameters are listed in Table 3. We adopt an $H$-band solar Vega magnitude $M_{\odot,H}=3.32$ from Table 1 of Blanton & Roweis (2007).

The gravitational potential of the galaxy is calculated from the summation of the MGE model components multiplied by a constant $M/L$ and the SMBH modelled as a point like mass, by following the equation in Cappellari et al. (2002). The circular velocity curve in the equatorial plane is then calculated from this by using the ; mge_circular_velocity procedure within the Jeans Axisymmetric Modelling (JAM) package444http://purl.org/cappellari/software (Cappellari 2008).

Our assumption of a spatially constant $M/L$ should be treated with caution given that the molecular gas disc (and potentially associated star formation) extends to the nucleus of NGC 3665 (Davis et al., 2014). The effect of any uniform star formation or young stellar population in the region of interest is implicitly subsumed into our adopted $M/L$, but a steep gradient could be problematic. However, neglecting centrally-concentrated star formation, and thus a decreasing $M/L$ with decreasing radius, effectively overestimates the stellar contribution to the total mass in the very centre. Our SMBH mass estimates are thus conservative. The effects of a potentially varying $M/L$ will be further explored in Davis & McDermid (2017).

3.2 Data Cube Model

Given the circular velocity curve obtained from the MGE formalism described above (Section 3.1), and the adopted molecular gas disc inclination (which is the same as that used to de-project our MGE models), we generate a model data cube using the Kinematic Molecular Simulation (KinMS) code of Davis et al. (2013)555https://github.com/TimothyADavis/KinMS. This assumes circular motions and a spatially uniform (but free) gas velocity dispersion. Instrumental effects such as beam-smearing and spatial and velocity binning are all taken into account by KinMS. For the properties of the molecular gas disc, we further assume that it has an exponential surface brightness profile with a void in the centre. The three free parameters describing the molecular gas distribution are thus its surface brightness scaling factor, radial scale radius and the radius of the central hole. The other free parameters required to fully describe a model (or, rather, to allow its comparison to real data) are the kinematic centre and position angle of the molecular gas disc in the plane of the sky and the galaxy systemic velocity. The total list of $11$ free parameters is given in Table 4.

3.3 Model Fitting

We use Bayesian analysis techniques to estimate the best-fit set of model parameters from our data cube, including the SMBH mass. Specifically, we utilize a Markov chain Monte Carlo (MCMC) method with Gibbs sampling to explore the parameter space. The number of iterations is set to $10^{6}$, and the first $5\times 10^{4}$ iterations are ignored as a burn-in phase. The method will be fully described in Davis et al. (2016), but we provide a short summary here. The aim is to obtain the posterior distribution of the $11$ model parameters: SMBH mass, stellar $M/L$, and the molecular gas disc kinematic centre, inclination, position angle, systemic velocity, velocity dispersion, integrated flux (CO surface brightness scaling factor), radial scalelength and void radius.

The region of the cube used for fitting covers the entire CO emitting region, and is defined to be $13\aas@@fstack{\prime\prime}0\times 13\aas@@fstack{\prime\prime}0$, the centre coinciding with the core position observed by Liuzzo et al. (2009). Velocity channels are from $-380$ to $380$ km s*-1* with respect to the systemic velocity of $2069$ km s*-1*.

We use a logarithmic likelihood function based on the $\chi^{2}$ distribution, calculated by comparing the CO distribution in each channel of the data cube with that in the model. As our data are approximately Nyquist sampled spatially, the synthesised beam induces strong correlations between neighbouring pixels in the data cube. The likelihood function we use, $\exp(-\chi^{2}/2)$, takes this into account by including the full covariance matrix when calculating the $\chi^{2}$. As the condition number of the covariance matrix itself is large, we do not invert it directly to calculate the likelihood, but instead introduce a modified Cholesky factorization step to avoid loss of numerical precision when calculating the inverse. The observational error on the flux in each pixel is set to the rms noise of the data cube evaluated in the central regions, in channels where no emission is detected.

We use flat priors for all the fitted parameters, within certain ranges. The prior distributions used and the posterior distributions returned are summarised in Table 4. The posterior distributions are also shown with greyscales in Figure 3. A comparison of the data and best-fit model moment 0 and moment 1 maps is shown in Figure 4. An analogous comparison of the channel maps is shown in Figure 5.

3.4 Model Results

We take the best fit of each model parameter directly from our Bayesian analysis, as the parameter value with the smallest $\chi^{2}$ in our multi-dimensional parameter space. For example, the SMBH mass is determined to be $5.75\times 10^{8}$ $M_{\odot}$, consistent with the predicted value from the known stellar velocity dispersion ($216$ km s*-1*; Cappellari et al. 2013) and the $M_{\rm BH}$–$\sigma$ relation of McConnell & Ma (2013).

The error bars of each model parameter are computed as percentiles of the posterior. Specifically, the $1\sigma$ confidence limits are the 15.9th and 84.1th percentiles of the posterior, while the 3sigma limits are the 0.14th and 99.87th percentiles of the posterior (see the grey-shaded regions in the likelihood contours of Figure 3).

Using this procedure, the SMBH mass and stellar $M/L$ are measured to be $5.75^{+1.49}_{-1.18}\times 10^{8}$ $M_{\odot}$ and $1.45\pm 0.04$ $(M/L)_{\odot,H}$, respectively, at the $3\sigma$ confidence level. The reduced chi-square ($\chi^{2}_{\rm red}$) for the best fit is $0.75$, indicating a good fit. The inclination angle is measured to be $69\degr\!\!.90\pm 0.61$ under the particular morphology of the molecular gas disc assumed. See Table 4 for the other best-fit parameters describing the molecular gas disc.

This SMBH mass gives an intrinsic SOI of $0\aas@@fstack{\prime\prime}3$, half of the synthesized beam. This SOI radius is slightly smaller than the radius of the cavity in the best fit, $0\aas@@fstack{\prime\prime}38$. The SMBH mass is thus constrained to $23\%$ at the $3\sigma$ level ($7\%$ at $1\sigma$), although we do not detect Keplerian motion in the very centre. Radial plot of the enclosed mass (Figure 6) shows that the best-fit SMBH mass is $6.98$ times the enclosed stellar mass within $0\aas@@fstack{\prime\prime}3$ (SMBH SOI), $3.73$ times the stellar mass within $0\aas@@fstack{\prime\prime}38$ (molecular void radius) and $1.71$ times the stellar mass within $0\aas@@fstack{\prime\prime}6$ (synthesized beam size). An uncertainty for stellar mass within the SOI is less than $2.0\times 10^{7}M_{\odot}$. We consider this to be negligible for the SMBH mass error. We further investigate our SMBH mass error budget in Section 4.3, focusing on possible systematic effects.

4 Discussion

4.1 CO Morphology

As NGC 3665 harbours an AGN, the molecular gas could be affected by its presence and/or that of the SMBH itself. With the highest angular resolution achieved so far, our CO observations have revealed that the molecular gas disc in NGC 3665 has an unresolved central hole, the position of which coincides with that of the radio core (detected with VLBI observations; Liuzzo et al. 2009) and the peak of the stellar surface brightness (identified from HST observations).

The absence of molecular gas within a radius of $\approx 0\aas@@fstack{\prime\prime}4$ or $\approx 65$ pc suggests that some mechanism may be dissociating the molecular gas or preventing it from forming or accumulating in the very centre of the galaxy. The dissociation of molecules generally has two main causes, AGN activity and UV radiation from young stars. While the former is clearly a possibility, the latter is unlikely to be significant in NGC 3665 as the star formation density is low (Davis et al., 2014). Dynamical effects can also affect the distribution and survival of molecular gas. Shocks and resonances are obvious possibilities, but it may also be that the strong shear expected near the SMBH (where the circular velocity curve varies with the radius $R$ as $R^{-1/2}$) can destroy molecular clouds, where most molecules are generally found. Such dynamical mechanisms will be analysed further in future works. Evidently, the lack of molecular gas in the central hole also suggests that there is no current cold gas accretion in that region. Having said that, observations targeting different molecules and/or transitions are necessary to prove that the central hole is truly devoid of gas, and to explore possible mechanisms for void creation.

4.2 AGN Properties

As mentioned above, given NGC 3665’s stellar velocity dispersion (Cappellari et al., 2013), the best-fit SMBH mass is in agreement with the latest $M_{\rm BH}$–$\sigma$ relations (e.g. McConnell & Ma, 2013). However, we can also investigate whether this SMBH mass is consistent with the known properties of the AGN. A radio jet was detected early on in NGC 3665 (Parma et al., 1986; Nyland et al., 2016), with the jet axis almost exactly perpendicular to the major axis of the central molecular gas disc. While AGN jets are generally thought to emerge perpendicularly to their accretion discs, it is also commonly accepted that there need not be a connection with the orientation of the large-scale disc of the galaxy. While the latter is the case for low-luminosity AGN such as Seyferts in late-type galaxies, the jet and large-scale molecular gas disc in NGC 3665 clearly have a connection. The kinematic position angle of the molecular gas disc is determined to be $26\degr$, roughly perpendicular to the position angle of the jet, $137\degr$ (determined by drawing a line from the northeast to the southwest blob seen at 5 GHz by Nyland et al. 2016).

The X-ray luminosity of NGC 3665 ($L_{\rm X}$) was estimated to be $10^{40.1}$ erg s*-1* in the $2$–$10$ keV energy range, extrapolated from the total energy within the Chandra energy range of $0.3$–$8$ keV by assuming a power-law spectrum $N(E)=\alpha\,E^{1.7}$ and a value of $\alpha$ derived from the same observations (Liu, 2011). Comparing this X-ray luminosity to the Eddington luminosity calculated from our best-fit SMBH mass ($L_{\rm Edd}=1.27\times 10^{38}\,M_{\rm BH}$ erg s*-1* $M_{\odot}^{-1}=10^{46.8}$ erg s*-1*), we obtain an Eddington ratio of $\log(L_{\rm X}/L_{\rm Edd})=-6.73$. This relatively low Eddington ratio suggests radiatively inefficient flows, including powerful outflows such as the radio jet observed in NGC 3665 (Merloni et al. 2003 and references therein). We therefore conclude that the AGN properties of NGC 3665 do not conflict with the derived SMBH mass nor with the observed molecular gas kinematics.

4.3 Other Error Sources on the SMBH Mass

Several high angular resolution observations of kinematics, aiming to measure SMBH masses, show Keplerian upturns in galaxy centres (see, e.g., Miyoshi et al., 1995; Barth et al., 2016b). Keplerian motion is expected in a potential dominated by a SMBH, where observations resolve the SMBH SOI and some emission arises from within the SOI. Spatially resolving the Keplerian region naturally allows to constrain the SMBH mass to high accuracy, as only a point mass yields such a behaviour, but our data do not show clear Keplerian motion. A possible reason for this is that CO emission is not present in the SMBH’s vicinity, as our model reveals a central cavity of radius $0\aas@@fstack{\prime\prime}38$ (see Section 3.4). Even if that were not the case, the synthesized beam ($0\aas@@fstack{\prime\prime}60\times 0\aas@@fstack{\prime\prime}56$; see Section 2.2) could have smeared out the information from within the SOI ($0\aas@@fstack{\prime\prime}3$; see Section 3.4).

Nevertheless, we measure a SMBH mass with an uncertainty of only $23\%$ (at $3\sigma$ confidence). The stellar $M/L$ and inclination, that have a direct influence on the SMBH mass, also have small uncertainties. Here we therefore investigate possible systematic errors on these parameters, and then discuss other possible effects that could increase the error budget.

A potential reason for the small error on the stellar $M/L$ is the rather large fitting area, set to $13\aas@@fstack{\prime\prime}0\times 13\aas@@fstack{\prime\prime}0$ so as to include all CO emission. A large fraction of this area is dominated by the stars rather than the SMBH. The uncertainty on the stellar $M/L$ (and in turn the SMBH mass and inclination) may thus decrease as the number of constraints (i.e. the area) increases, irrespective of whether the fit is good or not, simply because the model looses its freedom to vary ($\chi^{2}$ would otherwise increase unacceptably). We therefore narrow down the fitting area to $4\aas@@fstack{\prime\prime}0\times 4\aas@@fstack{\prime\prime}0$, and repeat the fit described in Section 3. The SMBH mass is then measured to be $M_{\rm BH}=(5.37^{+1.24}_{-1.10})\times 10^{8}$ $M_{\odot}$ and the stellar $M/L=1.49\pm 0.05$ $(M/L)_{\odot,H}$, with an inclination angle $i=69\degr\!\!.64\pm 0.77$, all at $3\sigma$ confidence level. This best-fit parameter set yields $\chi^{2}_{\rm red}=1.20$. Surprisingly, with fewer data, the constraints on the SMBH mass ($19\%$ error), $M/L$ ($3\%$ error) and inclination ($1\%$ error, all at $3\sigma$ confidence) are just as tight as the original result. Comparing these results with the ones from the original fit ($M_{\rm BH}=5.75^{+1.49}_{-1.18}\times 10^{8}$ $M_{\odot}$, $(M/L)_{H}=1.45\pm 0.04$ $(M/L)_{\odot,H}$ and $i=69\degr\!\!.90\pm 0.61$, all at $3\sigma$; see Section 3.4), we nevertheless notice a possible systematic error on the stellar $M/L$, but not on the SMBH mass or inclination.

The systematic error on the stellar $M/L$ could come from the MGE model (stellar luminosity profile), that has no associated errors. The stellar luminosity profile becomes degenerate with the stellar $M/L$ when calculating the circular velocity, and thus clearly affects the uncertainty on the stellar $M/L$. Our MGE model also fluctuates depending on, for example, how we define the region of dust attenuation (see Section 3.1 and Figure 2 for the MGE model fitting). We thus use two more MGE models with different definitions of the region affected by dust, and fit to an area of $13\aas@@fstack{\prime\prime}0\times 13\aas@@fstack{\prime\prime}0$. We first create a mask to cover the right half of the HST image, divided with a line as shown in the left-hand panel of Figure 7. The unmasked region within a $6\aas@@fstack{\prime\prime}0$ radius is used to create the new MGE. The resulting SMBH mass is $M_{\rm BH}=5.4\times 10^{8}$ $M_{\odot}$ with $(M/L)_{H}=1.41$ $(M/L)_{\odot,H}$, yielding $\chi^{2}_{\rm red}=0.73$. Second, we mask all the pixels with negative values in an unsharp-masked version of the HST image, created using a Gaussian of FWHM $0\aas@@fstack{\prime\prime}2$. The mask and the unsharp-masked image are shown in the right-hand panel of Figure 7. By using the MGE fit created from the unmasked pixels, the SMBH mass and stellar $M/L$ are measured to be $M_{\rm BH}=5.50\times 10^{8}$ $M_{\odot}$ and $(M/L)_{H}=1.24$ $(M/L)_{\odot,H}$, yielding again $\chi^{2}_{\rm red}=0.73$. The best-fit values of the stellar $M/L$ are however beyond the statistical $3\sigma$ error in both cases.

Our tests thus indicate that systematic uncertainties on the SMBH mass are $\approx 0.4\times 10^{8}$ $M_{\odot}$, likely of the order of the quoted statistical $3\sigma$ errors. The systematic uncertainties on the stellar $M/L$ are estimated to be $\approx 0.2$ $(M/L)_{\odot,H}$, larger than the quoted statistical $3\sigma$ errors.

Here we do not explicitly consider the effects of potential spatial variations of the stellar $M/L$ onto the SMBH mass error budget. This could however lead to non-negligible uncertainties, and we plan to test for this in a future work. The uncertainties on the stellar $M/L$ and SMBH mass would also increase by considering e.g., galaxy distance and HST photometric zero point uncertainties.

We further comment that the $\chi^{2}$ in our fit is calculated by considering only the rms noise level in the cube as the observational error. Including other sources of error would possibly increase the error on the fitting parameters.

4.4 Rotation Curve Fit

CARMA and ALMA observations just recently revealed that molecular gas kinematics enables SMBH mass measurements in nearby galaxies (e.g. NGC 4526, Davis et al. 2013; NGC 1097, Onishi et al. 2015; NGC 1332, Barth et al. 2016a). However, when comparing to models, the first two studies only used a rotation curve extracted from a position-velocity diagram (PVD) taken along the galaxy major axis. The procedure can be briefly summarised: (1) determine the kinematic major axis and draw a PVD along it; (2) estimate the mean velocity at each position along the PVD to determine the rotation curve; (3) compare these velocity measurements with analogous ones made from a model rotation curve extracted from a model data cube (assuming a set of model parameters) in an identical manner; (4) identify the best-fit model parameters using a $\chi^{2}$ analysis throughout parameter space.

In this sub-section, we thus re-derive the SMBH mass of NGC 3665 by fitting only the rotation curve extracted from our data, to verify that the value derived is in agreement with the full data cube fit of the previous sections. First, we fixed all molecular gas disc parameters to those obtained from the full data cube fit (see Table 4), leaving only the SMBH mass and stellar $M/L$ as free model parameters. We then extracted the kinematic major-axis PVD from our data cube, and fit a Gaussian to the line-of-sight velocity distribution at each position to determine the mean velocity at that position, only keeping for the model fit measurements with a signal-to-noise ratio $S/N>3$. The uncertainty at each position was set to the root mean square of the channel width of the observation and the FWHM of the fitted Gaussian. If the Gaussian FWHM was too large, generally indicating a bad fit, we also excluded that position from the model fit. A total of $55$ mean velocity measurements were thus made along the extracted PVD, creating an observed rotation curve along the kinematic major-axis of NGC 3665.

Model data cubes are generated for a range of SMBH mass and $M/L$ values, and model rotation curves extracted in the same manner as for the data. These are then compared to the data in a $\chi^{2}$ manner. Figure 8 shows $\chi^{2}$ contours in the parameter space. By taking the average of upper and lower edge of the lowest contour level ($\chi^{2}_{\rm min}+53$), we first determine the stellar $M/L=1.48\pm 0.01$ $(M/L)_{\odot,H}$ (grey shaded area in Figure 8) as a range to extract the $\chi^{2}$ distribution of the SMBH mass. The $\chi^{2}$ values are then marginalized over the $M/L$ range, and a polynomial fit to the $\chi^{2}$ is determined. $\Delta\chi^{2}\equiv\chi^{2}-\chi^{2}_{\rm min,fit}$ distribution of $M/L=1.48\pm 0.01$ $(M/L)_{\odot,H}$ (grey circles), marginalized values (black crosses) and the polynomial fit (black line) are shown in the left panel of Figure 9. The minimum $\chi^{2}$ for the polynomial fit along the SMBH mass axis, thus $\chi{2}_{\rm min,fit}$, is realized with the SMBH mass of $6.8\times 10^{8}M_{\odot}$. A polynomial fit to the $\chi^{2}$ distribution along the $M/L$ is then determined with the SMBH mass fixed to $6.8\times 10^{8}M_{\odot}$ (along the vertical black line in Figure 8), as shown in the right panel of Figure 9. The uncertainties on each parameter were estimated from the $99.7\%$ confidence interval (i.e. the parameter values with $\chi^{2}\leq\chi^{2}_{\rm min}+9$). The best-fit SMBH mass and stellar $M/L$ are then $6.8^{+1.1}_{-1.0}\times 10^{8}$ $M_{\odot}$ and $1.47\pm 0.02$ $(M/L)_{\odot,H}$, respectively, with a $\chi^{2}_{\rm min}$ of $129.8$ and $\chi^{2}_{\rm red}=2.45$ (the number of degrees of freedom here is $53$, allowing for the two free parameters). We note that the minimum $\chi^{2}_{\rm red}$ is above unity, suggesting the existence of a better fit with other molecular gas disc parameters. A scatter seen in the $\chi^{2}$ distributions (see Figure  9) can produce a systematic error for each parameters, but we do not further discuss about this here. We do not expect the systematic error of SMBH mass to be larger than $0.4\times 10^{8}M_{\odot}$ and $M/L$ to be larger than $0.2$ $(M/L)_{\odot,H}$, both of which is derived in Section 4.3. We stress that the SMBH SOI is not resolved in our observations, and the Keplarian motion is not seen in our PVD. Resolving the SMBH SOI can narrow down the error budget of the SMBH mass in this rotation curve analysis.

A comparison of the observed and best-fit model PVD and their associated rotation curves is also shown in Figure 10. The two parameters for each panel are $M_{\rm BH}=0$ with a best-fit $M/L=1.48$ $(M/L)_{\odot,H}$ (left panel; $\chi^{2}_{\rm red}=9.73$), the best-fit $M_{\rm BH}=6.8\times 10^{8}$ $M_{\odot}$ and $M/L=1.47$ $(M/L)_{\odot,H}$ (middle panel; $\chi^{2}_{\rm red}=2.44$) and an overweight SMBH with $M_{\rm BH}=3.4\times 10^{9}$ $M_{\odot}$ and $M/L=1.26$ $(M/L)_{\odot,H}$ (right panel; $\chi^{2}_{\rm red}=40.01$).

The SMBH mass and stellar $M/L$ derived from the rotation curve fit are both consistent with those obtained from the full data cube fit, despite the very small uncertainties on the $M/L$ ratios. In fact, the uncertainties derived are also very similar. The two fitting methods therefore appear equally useful a priori. This is probably because the rotation curve fit extracts the crucial part of the data cube, most closely related to the mass of the SMBH and stellar body. However, we recall that the rotation curve fit has only $2$ free parameters, compared to $11$ for the data cube fit. The other parameters, taken from the data cube fit, were fixed during the rotation curve fit. The data cube thus contains a lot of information useful to constrain many other model parameters associated with the structure of the molecular gas disc itself (see Section 3). In addition, a major disadvantage of the rotation curve fit is that some parameters can be degenerate with the SMBH mass or stellar $M/L$, and therefore must be independently constrained (e.g. the disc inclination, that directly affect the $M/L$ and to a lesser extent the SMBH mass). Extracting the appropriate PVD also requires a well-defined kinematic major axis, that is difficult to specify when the velocity field shows warps and/or kinematic twists. A middle ground is to fit an extracted velocity field. For example, Neumayer et al. (2007) considered warped and misaligned ellipses to model the ionised gas velocity field of NGC 5128 (Centaurus A). In NGC 3665, however, there is no significant evidence for non-circular motions. Models with warped structures may nevertheless be required in the future, as highly detailed gas distributions gradually become available through higher angular resolution observations.

5 Conclusions

We present CARMA 12CO$(J=2-1)$ observations of the early-type galaxy NGC 3665 with $0\aas@@fstack{\prime\prime}59$ resolution. These reveal a regularly rotating molecular gas disc in the equatorial plane of the galaxy, with an apparent void within a radius of $\approx 0\aas@@fstack{\prime\prime}4$ or $\approx 65$ pc, potentially created by the known AGN.

Fitting the entire observed data cube of NGC 3665 with a model with free SMBH mass, stellar $M/L$, and numerous parameters describing the structure of the molecular gas disc, we derive a SMBH mass of $M_{\rm BH}=5.75^{+1.49}_{-1.18}\times 10^{8}$ $M_{\odot}$ and a stellar $M/L$ of $(M/L)_{H}=1.45\pm 0.04$ $(M/L)_{\odot,H}$ at $3\sigma$ confidence levels (statistical error). The SOI of the SMBH is thus estimated to be $0\aas@@fstack{\prime\prime}3$, which is half of the synthesized beam. The central hole in the gas disc ($0\aas@@fstack{\prime\prime}38$) limits our gas tracer to radii outside the SOI, and all of the dynamical constraints on the SMBH mass come from model fitting to data outside of the SOI, where stellar mass is dominant. This SMBH mass is in agreement with that estimated from the latest $M_{\rm BH}-\sigma$ correlations, and appears consistent with the known AGN properties of NGC 3665, such as its radio jet and X-ray luminosity. Systematic uncertainties on the stellar $M/L$ are estimated to be $\approx 0.2$ $(M/L)_{\odot,H}$, by considering smaller fitting regions and several different MGE models. We estimate the systematic uncertainties on the SMBH mass to be $\approx 0.4\times 10^{8}$ $M_{\odot}$, which is within the statistical $3\sigma$ error. The full data cube fit also yields a SMBH mass consistent with that derived from a fit to the rotation curve only, but it opens the door to SMBH mass measurements in sources with significantly more complex molecular gas discs.

This work is only the fourth SMBH mass measurement using molecular gas kinematics, following measurements in two other lenticular galaxies and one barred spiral (Davis et al., 2013; Onishi et al., 2015; Barth et al., 2016a). This method has thus now proven its usefulness to derive SMBH masses in various types of galaxies. It offers exciting prospects to both calibrate SMBH masses measured with other methods, and to simply increase the number of galaxies with reliable SMBH masses. Further investigations comparing SMBH masses measured using other methods (stellar kinematics, ionised-gas kinematics and/or megamasers) will be required before a proper discussion of potential systematic differences between the different methods is possible.

Acknowledgements

The authors thank the referee, Aaron Barth, for constructive comments and useful suggestions. Part of this study was financially supported by JSPS KAKENHI grant number 26*368. KO acknowledges support from a Royal Astronomical Society. MB was supported by STFC consolidated grant ‘Astrophysics at Oxford’ ST/H002456/1 and ST/K00106X/1. MC acknowledges support from a Royal Society University Research Fellowship. TAD acknowledges support from a Science and Technology Facilities Council Ernest Rutherford Fellowship. Support for CARMA construction was derived from the states of California, Illinois and Maryland, the James S. McDonnell Foundation, the Gordon and Betty Moore Foundation, the Kenneth T. and Eileen L. Norris Foundation, the University of Chicago, the Associates of the California Institute of Technology and the National Science Foundation. CARMA science was supported by the National Science Foundation under a cooperative agreement, and by the CARMA partner universities. This research made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work is also based on observations made with the NASA/ESA HST, 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. Data analysis was partially carried out on the common use data analysis computer system at the Astronomy Data Center, ADC, of the National Astronomical Observatory of Japan.