RESOLVING THE DYNAMICAL MASS OF A z ∼ 1.3 QUASI-STELLAR OBJECT HOST GALAXY USING SINFONI AND LASER GUIDE STAR ASSISTED ADAPTIVE OPTICS

The tight correlations between the properties of supermassive black holes (BHs) and their host galaxies (Kormendy & Richstone 1995; Magorrian et al. 1998; Gebhardt et al. 2000; Ferrarese & Merrit 2000; Tremaine et al. 2002; McLure & Dunlop 2002; Marconi & Hunt 2003; Häring & Rix 2004; Graham 2004; Kormendy & Bender 2009; Gültekin et al. 2009; Jahnke et al. 2009; Merloni et al. 2010; Graham et al. 2011) are a powerful tool for probing their relative and possibly interdependent evolution and growth across cosmic time. Since the hierarchical assembly in conjunction with the central limit theorem can explain the existence of these correlations (Peng 2007; Hirschmann et al. 2010; Jahnke & Maccio 2011), evolution in the correlations may solely be due to the relative growth speeds of the BH mass and the corresponding total and bulge stellar masses of the host galaxy as a function of redshift. Hence, evolution in these correlations is a key diagnostic of the growth channels and fueling mechanisms of BHs, active galactic nuclei (AGNs), and their host galaxies.

Overall, the current picture is one of evolution in the various correlations clearly being strong at z > 2 (Walter et al. 2004; Peng et al. 2006a, 2006b; Somerville 2009; Decarli et al. 2010). A number of studies also find a strong evolution signal at low redshifts (Woo et al. 2006; Treu et al. 2007; Bennert et al. 2010), particularly in terms of correlations involving the bulge properties of the host galaxies: it is generally recognized that BHs of a given mass are hosted by increasingly less massive spheroids/bulges at larger values of z. However, at lower redshifts the majority of studies show a mild (by a factor ≲2–3) or absent evolution, particularly where the total stellar mass of the host galaxy is concerned. For example, for AGNs at 1 ≲ z ≲ 2 Merloni et al. (2010) observe an offset of a factor of two to three from the z = 0 scaling relation, and note that given the large scatter, the data could still be consistent with zero evolution at the 2σ level. And while Decarli et al. (2010) observe an overall evolution of up to a factor of seven between z = 0 and z = 3, the sources observed at intermediate (1.0 < z < 1.5) and low (z < 1) redshifts are consistent with modest (0.3 dex) and no evolution, respectively.

Beyond the nearby universe, observational estimates of BH and galaxy masses are subject to considerable uncertainties. Furthermore, different measures of galaxy mass (i.e., bulge mass, stellar mass, and total mass) are frequently used, necessitating additional care in the interpretation of results from different studies. The major uncertainties of previous work in this field affect both large and small samples equally, and arise from the accuracy with which the masses of BH and host galaxy (total and/or bulge mass) can be determined. Outside the local universe, BH masses can only realistically be obtained for currently accreting type 1 AGNs (i.e., quasars or Seyfert type 1 galaxies). For such objects, M_BH can be constrained from the combination of the UV-luminosity and the broad emission line widths from ionization models which assume virialized motion and are in turn calibrated via reverberation mapping on local AGNs with a statistical uncertainty within a factor of three to four for Mg ii and C iv (Vestergaard 2002; McLure & Jarvis 2002; Shields et al. 2003; Greene & Ho 2005; Vestergaard & Peterson 2006; Woo et al. 2006; Kollmeier et al. 2006; Salviander et al. 2007; Netzer & Trakhtenbrot 2007; Treu et al. 2007; McGill et al. 2008; Onken & Kollmeier 2008), or less for Balmer lines (Denney et al. 2009; Shen & Kelly 2010), plus any further systematic uncertainty associated with the calibration between reverberation mapping BH masses and dynamical BH masses at z = 0.

However, obtaining accurate (total and/or bulge) dynamical masses for the quasar host galaxy, the second key ingredient of the correlation in question, is very difficult for objects beyond z > 0.5. The presence of an active nucleus makes stellar velocity dispersion measurements from absorption lines infeasible, and the estimated velocity dispersions derived from narrow emission line widths (Shields et al. 2003, 2006; Salviander et al. 2007; Ho et al. 2008; Riechers et al. 2008, 2009) are fraught with uncertainty (Ho 2007; Greene et al. 2009). Inferring stellar masses from imaging-based measurements of the host galaxy stellar luminosity is another possible approach used by a number of previous studies (Borys et al. 2005; Peng et al. 2006a, 2006b; Jahnke et al. 2009; Somerville 2009; Merloni et al. 2010) but, as the mass-to-light ratios deduced from the UV/optical spectral energy distributions of star-forming galaxies are severely model-dependent, inferred stellar masses are non-unique/uncertain. At redshifts z > 4, practical mass determination is limited to methods such as the determination of dynamical masses via CO observations (e.g., Walter et al. 2004). Direct measurement of a host galaxy's dynamical mass is free of many of the problems associated with other methods and can equally be applied to observations at lower redshifts. Aside from millimeter observations of CO emission, dynamical masses can also be measured from interstellar medium (ISM) emission lines in the rest-frame optical.

At redshifts of z > 1, Hα is the most prevalent powerful emission line that can be used in practice as a kinematic tracer in AGN host galaxies. Even so, successfully probing the geometry and kinematics of the Hα-emitting gas requires integral field near-IR (near-IR) spectroscopic observations with a spatial resolution of not more than a few kiloparsec (i.e., ≲ 05). A sufficiently high (narrow-line) luminosity is needed, but the Hα luminosity of galaxies undergoing star formation at the average rate is easily sufficient—along with the presence of favorably orientated and ordered velocity structures (also likely for a significant proportion of objects). As a pilot study for direct dynamical host galaxy masses for z > 1 quasars, we present and model near-IR integral field spectroscopy of the z = 1.3118 quasar J090543.56+043347.3, using SINFONI (Bonnet et al. 2004; Eisenhauer et al. 2003). We intend to present the results of similar observations of our full sample in a subsequent paper; interim analysis suggests that the required flux, orientation, and ordered velocity field constraints are met for approximately two-thirds of the objects targeted. A standard cosmological model with Ω = 0.27, Λ = 0.73, and H₀ = 71 kms⁻¹Mpc⁻¹ is assumed throughout; this leads to a spatial scale of 8.45 kpc arcsec⁻¹ at this redshift.

Our target, the z = 1.31177 SDSS quasar J090543.56+043347.3 (Schneider et al. 2007; Shen et al. 2008), was selected on the basis of the presence of a nearby guide star (12th magnitude at a distance ∼43) to be used for tip–tilt correction. This object was observed using the near-IR integral field spectrometer SINFONI on the Very Large Telescope of the European Southern Observatory (ESO; Bonnet et al. 2004; Eisenhauer et al. 2003) and the PARSEC laser guide star adaptive optics, on the night of 2009 December 15, with an average airmass of ∼1.16 and integral field unit image seeing of 047 FWHM. Six pointings of 600 s were obtained, in two repeating patterns of target-sky-target in order to facilitate the removal of sky background emission, with a sky-frame declination offset of 30'' from the target coordinates. The H-band filter was used with a grating central wavelength of 1650 nm, giving a spectral coverage of 1.45–1.83 μm and a spectral resolution of 3000. The spatial scale of each of the 32 SINFONI slitlets was 100 mas; these are imaged onto 64 detector pixels apiece, giving a field of view of 3'' × 3'' for each pointing.

Analysis of SINFONI data requires a painstaking treatment of cosmic ray hits and bad pixels, neither of which can be adequately managed by the existing pipeline reduction routines alone. We use IRAF routines to identify bad pixels in each raw frame of data prior to applying the instrument pipeline, creating a mask file for each individual raw frame. Otherwise, the initial ∼2%–3% of bad pixels flagged as NaN in each raw frame would breed disproportionately as the raw data are straightened, stretched, shifted, and combined within the pipeline, resulting in an unrealistic fraction of the final data cube containing NaN values (typically ∼30%–50%). Ideally, we do not seek to remove signal from any pixel with some level of contamination by bad data, but rather to trace the propagation of bad pixels and thus gauge the reliability of each data element within the final data cube. From this point onward we work with two different data sets. First, we use our final bad pixel mask to linearly interpolate over all flagged bad pixels in the observed data. Although these flagged pixels will contribute flux to their neighbors, they will only do so with locally realistic count values. Hence, we do not invent data while preserving the integrity of the spaxel grid. Second, mask-frame images are generated and treated in the exact same manner as the observed data, thus allowing us to keep track of the behavior and growth of the bad pixels as the data cube is generated.

A further issue is the mislocation of slitlet edges by the ESO SINFONI data reduction pipeline. To counter this, we apply small pixel shifts to the supplied WAVE, MAP calibration images, recovering two to three of the 3–4 columns of data on each side of the field of view, i.e., 10% of the field, which would otherwise be lost. The outermost columns of each slitlet which cannot be reliably wavelength calibrated remain masked out. Each raw frame and its associated mask are then individually converted from a two-dimensional image into three-dimensional data cubes using the "sinfo-rec-jitter" routine in the pipeline and the improved skyline subtraction IDL routines of Davies (2007). Our bad pixel mask data cubes contain the proportional contribution of bad pixels to each spaxel in the cube. For spaxels where the proportional contribution of bad pixels is larger than 50%, that particular spaxel in the corresponding observational data cube is converted to a value of NaN; such spaxels do not contain reliable information. This accounts for a total of approximately 6% of spaxels in our final combined data cube. The remaining spaxels in the data cube are either completely uncontaminated by the presence of nearby bad pixels (true for ∼50%–70% of spaxels) or are dominated by emission from clean spaxels (∼25%–45%). As described earlier, in this latter case the spaxel flux is derived from an initial image where the flagged bad pixels/cosmic ray hits have been replaced by linear interpolation over the surrounding clean data. As these spaxels originate from a majority of good data pixels, with a minority contribution from proportionately scaled locally plausible flux values, they are judged reliable and left as they are. This done, the individual data cubes are straightened and aligned via a process of Gaussian centroiding, sub-pixel linear interpolation, and shifting of each wavelength slice.

Our final data cube is produced by mean combining the individual aligned cubes on a spaxel-by-spaxel basis, ignoring masked NaN values, and 3σ clipping any value which deviates significantly from the typical values found at that spaxel position and those of the neighboring spaxels across all individual aligned cubes. Using SINFONI observations of the standard stars Hip047903 (G1V) and Hip047235 (G3V), and template spectra (Pickles 1998) for their specific spectral types scaled to the appropriate magnitudes in the H band, we determine the relative flux response for each wavelength slice across the data cube and apply the derived correction (accurate to ∼12% in absolute flux values) to the cube data.

3.1. The Spatially Resolved Narrow Hα Emission

In order to reveal the underlying kinematics of the host galaxy, we need to extract the extended narrow-line emission from our data cube. We first remove the continuum emission in each spaxel via a least-squares fit to a second-order polynomial (excluding the wavelength region covered by the broad Hα emission) and then turn to the more difficult task of separating the line emission of the host galaxy from that of the AGN itself. As the broad lines originate only from the (unresolved) immediate vicinity of the AGN, the observed spatial variation of broad-line intensity in our data cube can be used to accurately constrain the point-spread function (PSF) for our data (as demonstrated by Jahnke et al. 2004). For each spaxel in our data cube, we isolate the wavelength ranges of the spectrum, which contain strong broad-line emission, but without contamination by skyline residuals, to either side of the spatially extended but spectrally compact narrow-line emission (Figure 1(a)). Using the emission on either side of broad Hα, we subtract any residual background continuum level and sum up the remaining flux in the selected regions of the broad line. This provides us with a relative measure of the broad-line flux on a spaxel-by-spaxel basis in the form of a PSF image which accurately traces the relative distribution of light from the quasar point source at all wavelengths within our data cube (Figure 1(b)). The quasar spectrum itself is well described by the spectrum of the peak spaxel in our data cube, which is dominated by AGN emission alone. We construct a quasar-emission data cube by scaling the quasar spectrum by the relative flux in each pixel of our PSF image. Subtraction of this data cube from our calibrated data cleanly and accurately removes both the broad and narrow Hα emission due to the AGN, and results in a new data cube which contains only host galaxy emission. Non-negativity within the wavelength region covered by the broad Hα emission is not enforced. For 20 spaxels close to the center of the PSF, errors in the subtraction of the highly luminous quasi-stellar object (QSO) emission dominate the resulting spectra and in a few cases result in mild oversubtraction of the QSO emission. We do not retrieve reliable information on the narrow Hα emission from these spaxels.

Zoom In Zoom Out Reset image size

For each spaxel of the QSO-subtracted data cube, we then fit the data with a linear continuum and a simultaneous combination of three Gaussians of the same width and fixed relative centroid but varying strength for the Hα and [N ii]6548 Å, 6584 Å emission lines. The spectra are converted to the source rest frame using the canonical Sloan Digital Sky Survey Data Release 3 (SDSS DR3) redshift of z = 1.31177. Figure 2 displays example fits to three different spaxel spectra in the QSO-subtracted data cube from different regions of the field of view. The results are then masked to exclude unreliable fits to noise rather than to genuine line emission, i.e., excluding the outermost rows and columns of data where no genuine emission features are observed, spaxels where the fit results in extreme velocities (−100 < v < 600 km s⁻¹ without subtraction of any further systemic velocity at the assumed redshift of 1.31177, confirmed as fits to noise features via visual inspection), line centroid uncertainties >40 km s⁻¹, and Hα fluxes of <10 counts or uncertainties of >3 counts per binned spaxel. Figure 3 displays the extracted Hα flux image (Figure 3(a)), line width image (Figure 3(b)), velocity field after correcting the systemic velocity (Figure 3(c)), and [N ii]6584 Å/Hα line ratio map (Figure 3(d)), overlaid with the PSF contours from Figure 1(b). A more in-depth analysis of the emission-line properties of this system is deferred to our study of a larger sample of similar objects. In our subsequent modeling of the velocity field, we further mask the outermost, low signal-to-noise regions, via boxcar-averaging the line flux in 9 × 9 pixel regions and excluding regions of low average signal. The remaining unmasked pixels are illustrated in the line flux image displayed in Figure 4(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Velocity Field Analysis

As is clear from Figure 3(c), this source displays a complex velocity field. A quick visual inspection suggests the central region can be described well by ordered disk-like rotation of the gas, with a reversal of the gas velocity at larger radii along at least one arm-like feature, and possibly a second outside our field of view to the north. With a pixel scale of ∼0.4 kpc pixel⁻¹, the southern arm extends to a projected distance of approximately 20 kpc, suggesting a tidal origin for this feature as the most likely scenario. Although there are no radio observations made specifically of this object, the fact that it is not a FIRST (Faint Images of the Radio Sky at Twenty cm; Becker et al. 1995; White et al. 1997) radio source suggests that no strong radio jet is present in this source, and thus that jet-driven shocks/interactions are unlikely to be a significant factor for the observed kinematics of this object.

Our aim is to use this observed gas velocity field as a tracer of the overall gravitational potential of the galaxy. The fact that disk-like rotation is present in the gas velocity field, as we will show in more detail in the following sections, is independent of the geometrical structure of the underlying stellar body and makes no statements about the stars being in either a disk-dominated or bulge-dominated structure.¹ Gas disks can be commonly observed in elliptical galaxies (Oosterloo et al. 2007, 2010; Emonts et al. 2008), and it is often the case that their kinematics do not trace those of the underlying stellar population. Indeed, the velocity structures may be highly decoupled, both in the case of active (e.g., HE1029−1408 at z = 0.086, which displays low stellar rotation but substantial gas rotation; Husemann et al. 2010) or inactive (e.g., NGC 2768, NGC 4314, NGC 4526, and NGC 4546; Sarzi et al. 2006) galaxies.

In order to properly characterize the rotation of the central regions and to extract a dynamical mass, it is important that we can also simultaneously account for the contrasting peculiar velocity structures dominating the gas velocity field at larger radii (>5 kpc). However, the signal to noise of the line emission in the individual spaxels where these two velocity structures overlap is not sufficiently high to separately resolve more than one line profile, nor can multi-spaxel extractions allow the fitting of multiple line profiles without unavoidable degeneracies in the relative line widths, centroids, and strengths. We therefore choose to use multi-component fits to the full velocity field.

3.2.1. Separation of Emission from the Central and Outer Tidal Components

In the process of producing a model velocity field, it is essential to consider the strength of the contribution of line flux at different velocities from adjacent pixels during convolution with the observational PSF. The observed velocity in a given position is not simply the convolution of the model velocity field with the observational PSF centered on that location: as the line flux varies across the observational field of view, the proportional contribution of different pixels to each specific value in the output model velocity field will also vary. In order to correctly weight each velocity component of our model during PSF convolution, we also need to include a model of the relative pixel-to-pixel line flux—for both the central rotating disk and also the additional tidal arm features.

In producing a model for the velocity field of J090543.56+043347.3, we start by examining the two major velocity structures separately, followed by more extensive multiple-component modeling. We obtain additional constraints on these velocity structures via the application of harmonic decomposition (cf. van de Ven & Fathi 2010; Fathi et al. 2005). This treatment divides a velocity field and associated signal-to-noise map into elliptical annuli of differing semi-major axes but identical orientations, inclinations, ellipticities, and centroids, and iteratively extracts the velocity profile as it appears within each annulus as a finite number of harmonic terms. Velocity and higher-order moments can then be readily characterized, and the nature of any deviations from non-circular motion assessed. The same tools can also be used to model a velocity field within a succession of elliptical annuli while also varying the parameters listed above for each annulus.

We use this technique to confirm the centroid of our velocity field, and find that the central regions are indeed well matched by ordered rotation, oriented at an angle of ∼10° east of north, and inclined to the line of sight at an angle of 35° ± 15°. The results of fitting the velocity field with a single simple inclined rotating disk, of fixed position angle and inclination angle, and a velocity which varies with radius, are displayed in Figures 4(a)–(d). These plots illustrate the velocity field, model, residuals, and radial rotation profile of the model, respectively (note that each data point in Figure 4(d) is derived from all spaxels which contribute to the model elliptical annulus at that radius). Velocity increases rapidly in the central regions to a maximum of 220  ±  50 km s⁻¹ and then declines rapidly beyond a radius of ∼12 pixels due to the influence of the counter-motions in the tidal features (Figure 4(d)).

We also carry out similar simplistic modeling of the velocity field targeting the outer regions. The velocity field due to the tidal arm features is oriented at an angle of approximately 140° offset from the central structures, varies relatively monotonically with radius, and dominates the velocity field beyond a projected radius of approximately 5 kpc (= 12 pixels; Figures 4(e)–(g)). Full harmonic decomposition along each set of the defined ellipses implies that neither the central rotating structure nor the outer tidal arms show any strong evidence for more complex motions within the regions they respectively dominate; the inclusion of any further velocity structures in our multiple-component modeling is therefore not required.

While harmonic decomposition provides a rapid means of investigating the different types of velocity structure contained with the velocity field of this galaxy, and determining better constraints on the system's inclination angle, we still need to model both structures simultaneously, particularly in order to determine the underlying rotational velocities in the overlap region. We model the light distribution of the tidal arm features by first tracing the observed line flux with a logarithmic spiral (Figure 5(a)). The average cross-sectional flux profile is determined and applied (Figure 5(b)), and subsequently convolved with the PSF (Figure 1(b)) to create a good match to the observed flux (Figure 5(c)). We then model the residual flux (Figure 5(d)) with an inclined exponential disk (oriented to match the central ordered rotation region of the velocity field), also PSF-convolved (Figure 5(e)), to provide a first approximation to the flux profile of the rotating disk. (Note that the residual flux after tidal-arm subtraction was better described as an exponential disk than as an n ⩾ 2 Sérsic profile, as expected for a rotating gas structure.) Figure 5(f) displays our combined relative flux model for the disk and tidal arm features, unconvolved with the observational PSF. While this does not perfectly match the observed flux profile of the narrow-line emission, it does provide a very good approximation to the relative pixel-to-pixel contributions to the velocity field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.2. Velocity Field Model and Dynamical Mass

The velocity profile of the arm features is modeled as a linear change in velocity with location along the arms, defined from the outer regions free of any influence from the central rotating disk, and extrapolated inward from there. For the central rotating disk, we assume for our model velocity field a Persic universal rotation curve (Persic et al. 1996) as a generally valid example of the rotation curve for a disk. This can be described as

$$\begin{eqnarray} V(x) &=& V(R_{{\rm opt}}) \left\lbrace\vphantom{\frac{x^2}{x^2 + 1.5^2 \left(\frac{L}{L_{\star }}\right)^{0.4}}} \left(0.72 + 0.44 \rm log \frac{{\it L}}{{\it L_{\star }}} \right) \frac{1.97x^{1.22}}{(x^2+0.78^2)^{1.43}}\right.\nonumber\\ &&\left.+ \,1.6\ {\rm exp} [-0.4({\it L/L}_{\star })]\frac{x^2}{x^2 + 1.5^2 \left(\frac{L}{L_{\star }}\right)^{0.4}} \right\rbrace \rm km\ s^{-1},\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where x = R/R_opt. We consider models with 0.1 < L/L_⋆ < 5.0, "optical" radii R_opt (roughly three times the disk scale length, and defined by Persic et al. (1996) to include ∼83% of the light from an exponential disk; in practice our models trace the velocity field out to twice the optical radius) ranging from 0.11 to 15 pixels (with a plate scale of ∼0.42 kpc pixel⁻¹ at this redshift). We also allow the inclination and orientation of the rotation axis to vary within our previously determined constraints from harmonic decomposition and from a prior model using a rotating disk alone without simultaneously accounting for the tidal arm structures. For each pixel, the contribution from the two-model components within a 21 × 21 pixel region equivalent to the size of the observational PSF are weighted according to their relative fluxes and convolved with the PSF.

Our best-fit model is determined using a Levenberg–Marquardt least-squares minimization technique (Markwardt 2009) and has L = 0.25 L_⋆, R_opt = 10.7^+3.0_{− 1.6}, and an inclination angle of 47° ± 12°. While the formal 1σ parameter errors for our best-fit model are 0.2 pixels for R_opt and 11 for the inclination angle, we note that reasonable fits to the data can also be obtained with inclination angles ranging from 35° to 59° before gross misfitting of the data becomes clear in the model-subtracted residuals. As the assumed inclination angle has the greatest impact on the derived dynamical mass, the full parameter errors given above are derived after applying these limiting values for the inclination angle, and account for the loosest possible constraints on disk inclination for this source.

Figure 6(a) displays the deprojected rotation curve of the central gas disk component as a function of radius in the rest frame for our best-fit model. We also display deprojected rotation curves for our upper and lower limit mass models as defined by the parameter bounds above. Derived from these curves using the gas motions as a tracer of the overall potential, we also plot the enclosed dynamical mass as a function of radius (Figure 6(b)). This latter plot assumes that v_rot = v_circ. In reality, non-circular gravitational motions of the gas will also be present (e.g., streaming and turbulence in the ISM) and the derived dynamical mass profile is, strictly speaking, a lower limit. However, as is clear from Figure 3(b), where the line emission originates solely from the centrally rotating regions—the eastern and western quadrants of the centrally rotating structure, which are not subject to possible contamination by the tidal features—the observed line widths are at their lowest and in fact are barely resolved above the instrumental resolution (∼70 km s⁻¹). The intrinsic line dispersion due to non-circular gravitational motion in the gas will be minimal, and thus v_rot ≈ v_circ is a good approximation. In Figure 7, we display the masked data, best-fit model, and model-subtracted residuals for the velocity field.

We can reliably trace the galaxy rotation curve out to a distance of between 10 and 15 pixels, equivalent to a galactic radius of 4.2–6.3 kpc. (Note that the 70% encircled energy radius for the PSF of our data cube is ∼02, equivalent to a projected distance of ∼1.7 kpc at the redshift of this source.) Ideally, we would determine a dynamical mass with reference to the scale length of the stellar light in the host galaxy, but no such information is available for this object. (While we have determined a best-fit measure of R_opt ∼ 10.7 pixels for the gaseous emission for this galaxy, this does not necessarily have any bearing on the scale length of the stellar light, and should not be used as a proxy in this case.) Within a radius of 5.25 ± 1.05 kpc (i.e., 10–15 pixel radius), the enclosed dynamical mass of the host galaxy implied by our best-fit model parameters (including errors) is 2.05^+1.68_{− 0.74} × 10¹¹ M_☉.

3.3. Star Formation in the Host Galaxy

In addition to constraining the dynamical mass, the narrow Hα emission can itself be used to estimate the star formation rate (SFR) in the host galaxy. We measure a total narrow Hα line luminosity of 1.22  ±  0.21 × 10⁴³ erg s⁻¹ from the unmasked regions of our data cube (including the 12% uncertainty in the absolute flux calibration). A further maximum of 5% of the total line flux may have been missed in the central 20 spaxels where errors due to the subtraction of QSO emission dominate, and from the small number of other masked pixels within the high signal-to-noise regions of narrow-line emission, which we add to the upper error bound. The origin of the Hα line emission is not necessarily solely due to the presence of young stars: AGN activity and shocks may also play an important role in the spatially extended ionization of the gas. While we cannot place our data on a standard Baldwin, Phillips, & Terlevich diagram (Baldwin et al. 1981; Veilleux & Osterbrock 1987) without the addition of a second-line ratio, sources with log([N ii]/Hα) < − 0.5 are predominantly observed to be star forming rather than AGN-photoionized sources (e.g., Groves et al. 2006). Figure 8 displays the unmasked regions of our data cube shaded according to whether the line ratio log([N ii]/Hα) is < − 0.5 or > − 0.5, i.e., whether or not stellar photoionization is likely to be the dominant ionization mechanism, or whether the observed line emission is better considered as part of the extended narrow-line region of the AGN. Star formation dominates in the tidal features, but is equally important in some of the more luminous regions of Hα emission toward the center of the galaxy. If the Hα line emission is due solely to star formation activity alone, the observed line flux is equivalent to an SFR of 96⁺²²_{− 17} M_☉ yr⁻¹ using the scaling of Kennicutt (1998). As the Hα emission can be produced by AGN photoionization as well as via photoionization by young stars, this value represents a robust upper limit on the overall star formation in this system. In Figure 9(a), we display this cumulative SFR as a function of radius. In Figure 9(b), we display the cumulative SFR as a function of radius considering only the data with line ratio values of log([N ii]/Hα) < − 0.5. In this case, we observe an overall SFR of 50 ± 10 M_☉ yr⁻¹, derived from Hα emission that is unlikely to have been photoionized by shocks or by the AGN, providing a useful lower limit on this quantity. The average integrated SFR per square kiloparsec (projected) is displayed for all data and also for the data with log([N ii]/Hα) < − 0.5 in Figures 9(c) and (d), respectively. These plots clearly illustrate the importance of star formation throughout this galaxy, and that it is not restricted solely to either the inner or outer regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.4. Broad-line Emission and the QSO Properties

Taking the literature data for J090543.56+043347.3, the mass of the central BH has been estimated as 9.02 ± 1.43 × 10⁸ M_☉, based on the width of its Mg ii emission line and the local continuum luminosity (Shen et al. 2008). Additionally, this source has a bolometric luminosity of 3.10 × 10³⁹ W (Shen et al. 2008), corresponding to a factor of 0.25 L_EDD for a BH of this mass.

We can also determine a BH mass based on measurements of the broad Hα emission line. Using the formalism of Greene & Ho (2005)

$$\begin{eqnarray} M_{{\rm BH}} &=& \big(2.0^{+0.4}_{-0.3}\big) \times 10^6 \left(\frac{L_{\rm H\alpha }}{10^{42}\, \rm erg\, s^{-1}} \right)^{0.55 \pm 0.02}\nonumber\\\ms &&\times\left(\frac{{\rm FWHM}_{\rm H\alpha }}{10^3\, \rm km\, s^{-1}} \right)^{2.06 \pm 0.06}\,M_{\odot }. \end{eqnarray} \tag{ 2 }$$

We derive a spatially integrated QSO spectrum by combining the spectrum of the peak spaxel of our flux-calibrated and continuum-subtracted data cube (dominated by QSO emission) with the PSF profile derived from the broad-line flux across the field of view. The broad Hα line has a measured FWHM of 3160 ± 170 km s⁻¹, and an overall luminosity of 109 ± 13 × 10⁴²erg s⁻¹, resulting in a BH mass estimate of 2.83^+1.93_{− 1.13} × 10⁸ M_☉. We note that the Mg ii- and Hα-based BH masses differ by a factor of ∼3, and that the observed bolometric luminosity of Shen et al. (2008) is equivalent to a factor of 0.79L_Edd if we assume the Hα-derived BH mass, which is likely to be subject to a slightly lower overall systematic uncertainty.

The driving science question behind this work is to investigate the BH mass versus host galaxy mass relation at high redshifts, in as direct and unbiased manner as possible. Obtaining dynamical quasar host galaxy masses via narrow-line emission kinematics is an ideal way forward. Although the existence of complex velocity structures can potentially present difficulties for such a method, our pilot object J090543.56+043347.3 despite showing a non-simple velocity field can nonetheless be successfully modeled. Having obtained a measure of both BH and host galaxy dynamical mass for J090543.56+043347.3, in Figure 10 we contrast it with the canonical low-redshift M_BH versus M_{Bulge, Dyn} relation of Häring & Rix (2004) and also the more recent M_BH versus M_Stellar and M_BH versus M_{Bulge, Dyn} relations of Sani et al. (2011), based on Spitzer/IRAC data. The dynamical bulge masses of Häring & Rix (2004) plotted in Figure 10(a) are either derived via first confirming that the light profiles of the galaxies are bulge-dominated and then solving the spherical Jeans equation, or taken directly from the literature. This sample uses reliable preexisting BH mass estimates from the literature, derived from gas and stellar kinematics and primarily sourced from Tremaine et al. (2002). For the points on Figure 10(b), Sani et al. (2011) carry out two-dimensional image decomposition of the 3.6 μm data for their sample and then derive virial dynamical galaxy masses as M_{Bulge, Dyn} = 5σ²R_e/G. Once again, the BH masses for this sample are reliable estimates from the literature based on stellar or gas dynamics, or masers. Note that there is no significant difference between the regression lines for M_BH versus M_{Bulge, Dyn} and M_BH versus M_Stellar derived by Sani et al. (2011) from the 3.6 μ mass-to-light ratio including color-correction terms. Two separate points are plotted for the different Mg ii- and Hα-based BH masses of J090543.56+043347.3. We find that, using the Mg ii-based BH mass, the position of this object on either plot is a factor of ∼2 from the mean z = 0 relation, but well within the scatter. Using our Hα-based BH mass measurement, the position of this object lies directly upon the mean z = 0 relation. This no-evolution result is also consistent with the theoretical modeling of Jahnke & Maccio (2011), in which the correlation can be fully explained to have emerged at high redshifts from hierarchical assembly and merging statistics. Due to the drop both in merger rate as well as star formation and BH accretion rate densities in the universe since z ∼ 2, their explanation predicts the overall scaling relations to be largely in place by z = 1.

Zoom In Zoom Out Reset image size

It should be noted that in Figure 10 we compare the dynamical mass derived for a z = 1.3 galaxy with an unknown stellar morphology and a rotating gas disk with stellar masses and bulge dynamical masses in the local universe. Given this, it is therefore worthwhile to consider the likely subsequent evolution of this system over the intervening 8–9 Gyr. The major mechanisms for the growth and evolution of the BH and its host galaxy are BH accretion, star formation, merger activity, and disk-to-bulge reprocessing, which we now consider in turn.

The bolometric quasar luminosity of J090543.56+043347.3 (Shen et al. 2008) is equivalent to a BH accretion rate of ∼5.4 M_☉ yr⁻¹ (assuming mass is converted to luminosity with an efficiency of 10%). Quasar lifetimes are typically between 10⁷ and 10⁸ years (Yu & Tremaine 2002; Marconi et al. 2004; Hopkins & Hernquist 2009; Kelly et al. 2010). Taking the average of the two BH mass measurements for J090543.56+043347.3 (∼6 × 10⁸ M_☉) and a timescale of 5 × 10⁷ years, an accretion rate of ∼5.4 M_☉ yr⁻¹ would result in the central BH growing by ∼45% in total over the current active cycle of AGN-mode accretion. The absolute maximum possible growth factor for this BH (assuming (1) that it is observed at the very start of its current activity cycle, (2) a maximum accretion lifetime of 10⁸ years, and (3) the lower of the two BH mass estimates) is a factor of ⩽2, though we do stress that growth approaching this limiting value is unlikely. At the same time, it is statistically not very likely that this BH would undergo any further episodes of AGN activity between z = 1.3 and z = 0 (given the expected AGN fraction as a function of redshift; Bluck et al. 2011).

Considering star formation in the host galaxy, our observed Hα flux is consistent with an SFR of up to 100 M_☉ yr⁻¹ (Section 3.3). Only accounting for the line emission which is likely due to ionization by young stars rather than the AGN or shocks (particularly true of the tidal arms and the areas of high emission line intensity within the host galaxy: see Figure 3(d), Figure 9) gives an SFR of 40–60 M_☉ yr⁻¹. The observed SFR is significantly higher than average (e.g., Karim et al. 2011) for a galaxy at this redshift with a stellar mass approximating our derived dynamical mass (2.05^+1.68_{− 0.74} × 10¹¹ M_☉ within a radius of 5.25  ±  1.05 kpc). But this may not be surprising given the clear presence of highly star-forming tidal features (Figure 3).

As an aside, we note that triggering of the quasar activity via interactions/mergers is a possibility for this particular AGN, given the clearly discernable merger signature in the gas. The high projected SFR per square kiloparsec toward the center of the host galaxy (within a radius of 4 kpc; see Figures 9(c) and (d)) might also be suggestive of links between mergers, AGN activity, and circumnuclear star formation in the case of this source, or of a nuclear star formation ring such as that observed in the case of NGC 1097 (van de Ven & Fathi 2010). However, the triggering of AGN activity via mergers remains unproven and we are far from being able to, e.g., determine and compare the age of the AGN² with that of stars newly formed by a merger. Overall, there is very little direct observational evidence that galaxy mergers trigger AGNs even in individual cases at lower redshifts (e.g., Canalizo & Stockton 2000), while on a population basis major merging can clearly be ruled out as the dominating mechanism at z < 1 (Cisternas et al. 2011).

Assuming that the local scaling relations for BH mass versus galaxy stellar mass were already approximately in place at z ∼ 1.3, we note that mergers involving galaxies with similar M_BH/M_Total ratios will (statistically) lead to the same proportional growth for both BH and host galaxy, leaving the relation unchanged. It is expected that a stereotypical galaxy with the same mass as that determined for J090543.56+043347.3 would undergo at most one more major merger event (Robaina et al. 2010), and quite likely be restricted to more minor accretion events. Therefore, both the host galaxy and its central BH are likely to grow by a maximum of approximately 50% in mass by the present day due to merger activity alone (see, e.g., Hopkins et al. 2010; van Dokkum et al. 2010). Similarly, under the assumption that the observed scaling relations for M_BH versus M_Host are already in place at z = 1.3, merger activity alone would shift the exact placement of this object on Figure 10 in a direction roughly parallel to the observed mean z = 0 scaling relation.

Even though the gas and stars trace the same gravitational potential, the gas and stellar kinematics can be expected to differ. While the observed gas kinematics for the bulk of this system are relatively settled and consistent with a rotating disk, the stellar morphology cannot be inferred from the gas morphology, and we emphasize that it remains unknown. However, given the observational evidence for ongoing tidal activity, even if this galaxy is restricted in the future solely to minor accretion rather than major merger events, the observed dynamical mass at z = 1.3 can be readily reprocessed into a stellar bulge. This object is almost certainly the progenitor of a z = 0 bulge-dominated galaxy.

What bearing does the combination of all our observations have on the likely subsequent evolution of this system and on its position relative to the local scaling relations?

On the basis of the observations of J090543.56+043347.3, we can quantify the relative rate at which both the mass of the central BH and of the host galaxy's stellar population are growing at the epoch at which we observe them. The annual increase in BH mass of 5.4 M_☉ yr⁻¹, as derived from the quasar's bolometric luminosity, is just 10 times less than the observed (conservative) SFR of ∼50 M_☉ yr⁻¹. By contrast, the scaling relations predict a BH which is of the order ∼500 times less massive than its host galaxy. While the central BH is actively accreting, the host galaxy is forming new stars at a much more moderate rate relative to its overall stellar mass.

Of course, we know that the duty cycle of star formation for z ∼ 1.5 galaxies is far higher than the duty cycle of active BH accretion. Over a single AGN activity cycle lasting typically 10⁷–10⁸ years, the proportional increase in mass of the BH is expected to be ∼45%, as noted above. A similar proportional increase in the stellar mass of the host galaxy can be achieved if the observed level of star formation persists over a substantially longer timescale than the AGN activity: a galaxy of 2 × 10¹¹ M_☉ with an SFR of ∼50 M_☉ yr⁻¹ would increase in stellar mass by 45% over the course of ∼2 Gyr. In such circumstances, the BH mass versus host galaxy mass ratio of this object would remain consistent with the local relation. However, we note that the required star formation timescale of several Gyr is well beyond the typical merger timescales.

An alternate approach to constraining the expected future increase in stellar mass of J090543.56+043347.3 with time until z = 0 is to consider SFRs for the galaxy population as a whole. More specifically, following the cosmic average (e.g., Karim et al. 2011) the expected levels of star formation activity would result in a growth of typically a factor of about 1.5–3 in stellar mass between z = 1.3 and the present day, depending on the exact details of the subsequent star formation activity within the galaxy.

We can combine this empirical estimate of the stellar mass growth through star formation activity with the increase in mass due purely to the expected merger activity (∼50%; Hopkins et al. 2010). Overall, the host galaxy is very likely to more than double in mass (increasing by ≳ 0.7 dex in log(M_total)) between z = 1.3 and z = 0. As the majority of any disk stellar mass will have been reprocessed into the galactic bulge via galaxy mergers or tidal interactions, we would expect log (M_{total, z = 0}) ≈ log (M_{Bulge, z = 0}).

Similarly, we can combine the expected growth in BH mass due to AGN accretion and merger activity to predict the overall increase in M_BH by z = 0. These processes are likely to result in an overall increase of ≲ 0.25 dex in log(M_BH) for a system such as this; taking our less likely absolute upper limit for BH accretion, the maximum increase in M_BH is ≲ 0.7 dex.

Any change in position of this object by z = 0 relative to the local relation will therefore most likely be dominated by changes in the host galaxy mass. Overall, regardless of the exact details of the further evolution of this object, it is most likely that the average of the two points plotted for this object on Figure 10 will move to a position which is even closer to the best-fit z = 0 relations displayed in Figure 10, and in any case will remain consistent with the observed scatter.

In conclusion, we return to the observational aspects of the work presented here. We believe that we have presented a study that is ground-breaking in the respect that it demonstrated how adaptive-optics-assisted near-IR integral field spectroscopy can be successfully applied and used to determine the dynamical mass profile for the host galaxy of a typical quasar. The observed BH masses and overall dynamical mass of the host galaxy of this object at z = 1.3 are consistent with the observed present-day relations for M_BH versus M_{Bulge, Dyn} and M_BH versus M_Stellar, and aside from a potential conversion of disk to bulge mass through minor mergers or the likely ongoing major merger, no significant evolution in the BH versus stellar mass relationships between z = 1.3 and z = 0 is required or expected. Provided that one takes into account the minimum extended emission line flux requirements and the likelihood of favorable gas kinematics and orientation in the observed velocity field, this particular diagnostic route is widely applicable to other quasars in the redshift range 1 < z < 3 at the peak epoch of quasar activity. Analysis of a further seven objects in our pilot sample is currently underway. Although J090543.56+043347.3 itself displays a relatively impressive SFR, this diagnostic method is clearly feasible at much lower narrow Hα fluxes, while alternatives such as CO-based dynamical masses are most successful for extreme, atypical objects.

K.J.I. and K.J. are funded through the Emmy Noether Programme of the German Science Foundation (DFG). We thank Ric Davies, Chien Peng, and Arjen van de Wel for very useful discussions, Eleonora Sani for providing the z ∼ 0 data which appears in Figure 10(b), and the anonymous referee for detailed consideration of the manuscript.

Facility: VLT:Yepun - Very Large Telescope (Yepun)