DO THE MOST MASSIVE BLACK HOLES AT z = 2 GROW VIA MAJOR MERGERS?

As discussed in Section 1, a well defined selection function is necessary to understand the statistical biases present in any AGN sample. Previous contiguous-field morphological studies with HST (e.g., Grogin et al. 2005; Gabor et al. 2009; Cisternas et al. 2011; Kocevski et al. 2012) are substantially volume limited and thus lacking luminous quasars at the highest black hole masses (${M}_{\mathrm{BH}}\gt {10}^{9}\,{M}_{\odot }$). We are thus interested in extending these studies to test whether black hole accretion is merger driven in the highest-mass black holes at the highest redshift where HST studies can reliably diagnose merger signatures—z ≃ 2, where the rest-frame optical emission shifts into the WFC3 IR $F160W$ filter.

We selected quasars from the Sloan Digital Sky Survey (SDSS) 5th Data Release Quasar Catalog (Schneider et al. 2007), using virial black hole masses calculated from the Mg ii line by Shen et al. (2008). We required that the quasars be targeted using the uniform color selection of SDSS (TARGET_QSO_CAP or TARGET_QSO_SKIRT, see Richards et al. 2002, 2006), have a redshift in the range z = 1.9–2.1, and SMBH mass ${M}_{\mathrm{BH}}={10}^{9.3}\mbox{--}{10}^{9.7}\,{M}_{\odot }$. In order to ensure an unbiased sample of all optically luminous massive AGN, the quasars were selected blindly from the parent sample, with no further criteria based on spectral features, broadband colors, or detections at other wavelengths. Of this parent sample, we submitted 115 randomly drawn targets for an HST Cycle 19 SNAPshot survey (Program SNAP 12613, PI: Jahnke). Between 2011 October and 2012 September, 19 of these quasars were observed with the WFC3 infrared channel in the $F160W$ filter (rest-frame V-band). After the HST observations were completed, we examined existing data from the SDSS spectra and wide-area surveys. These data are summarized in Table 1 and discussed below.

Table 1.  Properties of Observed Quasars

Quasar	Redshift	L3000	FWHM	MBH	$L/{L}_{\mathrm{Edd}}$	${m}_{F160W}$	MV	Notes
(SDSS J)	z	$\mathrm{log}({L}_{\odot })$	Mg ii (km s−1)	$\mathrm{log}({M}_{\odot })$		mag	mag
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
081518.99+103711.5	2.021	12.72	3600	9.3	0.49	18.55	−27.49	⋯
082510.09+031801.4	2.035	12.91	4900	9.7	0.28	17.81	−28.27	⋯
085117.41+301838.7	1.917	12.54	3000	9.0	0.65	18.85	−27.08	⋯
094737.70+110843.3	1.905	12.66	6200	9.7	0.11	18.91	−27.01	C IV BAL
102719.13+584114.3	2.020	12.67	3700	9.3	0.57	18.67	−27.38	⋯
113820.35+565652.8	1.917	12.76	4300	9.5	0.30	18.12	−27.81	⋯
120305.42+481313.1	1.988	12.79	4100	9.4	0.27	18.18	−27.84	⋯
123011.84+401442.9	2.049	13.17	5400	9.9	0.28	17.41	−28.67	⋯
124949.65+593216.9	2.052	12.88	4100	9.5	0.35	18.11	−27.99	⋯
131501.14+533314.1	1.921	12.70	8400	10.0	0.10	18.45	−27.47	⋯
131535.42+253643.9	1.926	12.57	3200	9.1	0.71	18.50	−27.44	⋯
135851.73+540805.3	2.066	12.80	3700	9.4	0.54	18.27	−27.83	⋯
143645.80+633637.9	2.066	13.22	6800	10.1	0.22	17.06	−29.04	Radio-loud
145645.53+110142.6	2.017	12.79	3200	9.2	0.65	18.20	−27.86	⋯
155447.85+194502.7	2.091	12.71	5900	9.7	0.14	18.77	−27.37	⋯
215006.72+120620.6	1.993	12.66	5100	9.5	0.29	18.55	−27.46	⋯
215954.45−002150.1	1.963	13.33	4000	9.8	0.74	16.87	−29.12	Hot dust-poor
220811.62−083235.1	1.923	12.77	3300	9.2	0.63	18.50	−27.43	⋯
232300.06+151002.4	1.989	12.81	5400	9.7	0.19	18.13	−27.89	⋯

Note. Properties of z = 2 quasars. Columns 1–5 are from the catalog of Shen et al. (2011). Column 1: SDSS name, including the full sexagesimal celestial coordinates; Column 2: systematic redshift; Column 3: rest-frame 3000 Å luminosity; Column 4: FWHM of the broad component of the Mg ii line emission (median fractional error is 11%); Column 5: Mg ii virial BH mass (calibration of Shen et al. 2011); Column 6: Eddington ratio, derived using the Mg ii BH masses and the catalog bolometric luminosity; Column 7: $F160W$ observed magnitude from our HST observations (photometric error is $\simeq 0.05$ mag); Column 8: V-band absolute magnitude (k-corrected using the Vanden Berk et al. (2001) median quasar spectrum); Column 9: special notes or features.

Download table as:  ASCIITypeset image

Shen et al. (2011) classified the spectral properties of all the SDSS Data Release 7 quasars (Schneider et al. 2010), including the 19 in our HST program. In Table 1 we reproduce those data relevant to deriving the Mg ii virial BH masses—rest-frame luminosity at 3000 Å and FWHM of the broad Mg ii component. We note that Shen et al. (2011) also publish a revised calibration for Mg ii virial BH masses, different from the McLure & Dunlop (2004) calibration used by Shen et al. (2008) that went into our mass selection. The new calibration simply increases the scatter in BH masses—the sample is still representative of luminous AGN with the high BH masses and the median mass, ${M}_{\mathrm{BH}}={10}^{9.5}\,{M}_{\odot }$, remains unchanged. Table 1 also includes an estimate of the Eddington ratio for each quasar, calculated from the Shen et al. (2011) bolometric luminosity and Mg ii BH mass.⁸ One object was flagged as a C iv broad absorption line (BAL) quasar. We manually examined the SDSS spectra of all 19 objects and confirmed that it is the only BAL quasar in the sample. Eighteen of the quasars were covered by the FIRST 1.4 GHz radio survey (White et al. 1997, version 14Dec17) with a detection limit sufficient to assess radio-loudness—one object is radio-loud with ${S}_{1.4\mathrm{GHz}}=900$ mJy and R = 1670, and is core-dominated under the classification scheme of Jiang et al. (2007). In addition, one object is classified as hot dust-poor, defined as having a rest-frame mid-IR flux deficiency (Jiang et al. 2010; Jun & Im 2013), with a $2.3\,$^–$0.51\,\mu {\rm{m}}$ luminosity ratio of $\mathrm{log}({L}_{2.3}\,/{L}_{0.51})=-0.71$. One object out of 19 corresponds to an inferred fraction for each of these special object classes (BAL, radio-loud, hot dust-poor) in the range $0.035\mbox{--}0.15$ (beta distribution, 68% confidence interval, see Section 5.2 below). This is consistent with these fractions for the z = 2 luminous quasar population as a whole ($\simeq 10 \% \mbox{--}15$% BALs Gibson et al. 2009; Allen et al. 2011, 10%–15% radio-loud Jiang et al. 2007, $\simeq 2$% dust-poor Jun & Im 2013). Table 1 also includes the measured $F160W$ magnitude of each quasar from our HST imaging (see next section), and the rest-frame V-band absolute magnitude, calculated using the Vanden Berk et al. (2001) median quasar spectrum for the k-correction.

All 19 quasars were observed with the WFC3 infrared channel using the $F160W$ filter (broad H-band). None of the quasars had any existing HST imaging at rest-frame optical wavelengths. As a SNAP program, the observations of necessity have short integration times of 1597 s per target ($\lt 1$ orbit). Each observation was split into four exposures, dithered using the standard four-point half-pixel box pattern to improve the PSF sampling, and to assist in the rejection of bad pixels and cosmic rays.

Our data processing began with individual flat-fielded, flux-calibrated exposures delivered by the HST archive. The four exposures for each pointing were combined using the ASTRODRIZZLE software package (Koekemoer et al. 2002, 2013) with an output plate scale of 0060 per pixel and a pixfrac parameter of 0.8. For $F160W$ observations, this samples the PSF with 2 pixels per FWHM and provides relatively uniform weighting of the individual pixels. We used "ERR" (minimum variance estimator) weighting for the final image combination step. We also generated variance maps that include all sources of noise, including uncertainty from the CALWF3 count rate determination, by copying the WFC3 "ERR" arrays into the standard image arrays, and re-running the drizzle process using the same parameters and weighting scheme. The variance maps are a requirement for our analysis, since the count rate from the quasar point source is significantly higher than that of the sky, and underestimated errors can lead to problems with multi-component fitting, as described in Section 3.2 below.

Despite the short exposure time, the excellent sensitivity of WFC3 IR and low on-orbit near-infrared sky allow the images to reach $1\sigma$ limiting surface brightness in the range ≃24.0–24.5 mag arcsec⁻². This is sufficient to detect tidal features of major mergers between luminous galaxies, as illustrated by the WFC3 imaging of the most luminous, distorted z ≃ 2 merger in the GOODS/ERS2 field (Van Dokkum & Brammer 2010; Ferreras et al. 2012). We note that the true achieved surface brightness sensitivity is a function of distance from the quasar point source (due to shot noise from the removed point source). The sensitivity to low surface brightness features is discussed in detail in Appendix B.

As mentioned in Section 1, particular care must be taken in selecting quasar samples, since there is no unbiased selection method that captures the entire population. Furthermore, selecting non-random subsamples from even large survey catalogs may introduce additional biases, and observational constraints may bias sensitivity to features of interest. For our study, two kinds of bias are the most salient. First, if a universal evolutionary sequence such as that proposed by Sanders et al. (1988) exists, certain quasar selection methods may bias toward a particular phase within that evolution. Second, since we are interested in assessing evidence for mergers via morphological signatures, observational biases affecting sensitivity to those signatures are relevant, in particular surface brightness sensitivity and the rest-frame emitted wavelength. We discuss the observational biases of this program in detail in Appendix B.

Selecting optically luminous broad-line sources (Type 1 quasars) by definition selects objects where the central accretion disk is essentially unobscured along the line of sight. Highly obscured or very red (Type 2) sources missed by this selection may include objects similar to the unobscured sources but viewed from different angles, as well as fundamentally different objects that contain systematically more dust or different dust geometries. For luminous quasars such as those in this study (${L}_{\mathrm{Bol}}\gtrsim {10}^{47}$ erg s⁻¹), the obscured fraction due to anisotropy of the dust torus is estimated from Type 1 quasar spectral energy distribution (SED) modeling to be $\simeq 20 \% \mbox{--}50$% (Lusso et al. 2013). Mid-infrared quasar selection methods directly estimate the total luminous Type 2 fraction—i.e., the above torus-obscured quasars, but also those obscured due to fundamentally different dust geometries—to be $\lesssim 50$% (Donley et al. 2008; Assef et al. 2015). We thus expect $\lesssim 30$% of luminous quasars to be hosted in systems with much larger dust covering fractions than our sample, whether these represent the putative "buried" evolution phase with the quasar completely enshrouded in dust, or simply a sub-population of luminous quasars hosted in ULIRG-like dust-rich galaxies. This also sets an upper limit on a buried phase duration of at most 30% of the average quasar lifetime. We argue in Section 6.2 below that this timescale is insufficient for dynamical effects to erase the major merger signatures we are searching for.

Besides the implicit luminosity and unobscured line-of-sight constraints, there are potential biases from explicitly selecting quasars with high-mass black holes. In particular, since the (active) black hole mass function drops sharply at high masses, there could be concern that such high-mass quasars are preferentially near the end of their accretion phases. Although the volume density of active black holes drops at high mass, the Eddington ratio distribution function does not depend strongly on black hole mass (Schulze et al. 2015), so even high-mass black holes like those in our sample seem to have accretion episodes statistically similar to lower-mass black holes. We discuss merger signatures as a function of black hole mass in Section 6.2, but note here that even if the quasar duty cycle is near unity and these quasars grow continuously near the sample average $\langle L/{L}_{\mathrm{Edd}}\rangle =0.40$, 15 of the 19 could continue growing for $\gt {10}^{8}$ years before they exceeded the maximum mass of our selection region, on the long side for quasar lifetime estimates (Martini & Weinberg 2001; Yu & Tremaine 2002; Shen et al. 2007). Thus, although these high-mass black holes might be experiencing their last major growth phase, there is no a priori reason to assume we are observing the quasars late in the current active phase.

As a SNAP program, our HST data did not have dedicated observations of stars to measure the instrument and telescope PSF. The focus of all HST instruments is affected by the telescope's thermal environment, with changes in solar illumination resulting in de-space of the secondary mirror, the so-called "spacecraft breathing" effect (Bély et al. 1993; Hershey 1998). For applications requiring precise PSF matching, the exposure focal history is best matched by extracting PSF stars from the same images. However, the WFC3 infrared channel PSF has other aberrations (e.g., coma, astigmatism) that vary with position within the WFC3 field of view. Simulated WFC3 PSF models are currently poor matches to observations (e.g., Mechtley et al. 2012; Biretta 2012), so we chose to build a library of empirical PSF models from WFC3 archival data of high signal-to-noise ratio (S/N) stars near the center of the field of view, where all the quasar targets were observed.

To find similar exposures from which to extract PSF models, we searched the HST archive for all single-orbit $F160W$ observations using a 4-exposure dither pattern. We then identified all point sources falling within 05 of the WFC3 field of view center, and excluded known quasars or radio sources in the NASA/IPAC Extragalactic Database or were from HST programs specifically targeting AGN. We visually inspected the remaining PSF stars, and excluded any that were contaminated by background galaxies or whose flux distribution was significantly elliptical, i.e., were likely stellar binaries. This left us with a library of eight stars⁹ with high S/N, but which still had an accurate count rate determination in their cores. We created drizzled images and variance maps of these stars using the same plate scale and weighting scheme described in Section 2.2 above.

While a posteriori estimation of HST focus is possible from on-orbit thermal measurements using the HST Focus Model (Cox & Niemi 2011), there is yet another (essentially) degenerate source of PSF mismatch, namely the object SED through the broadband $F160W$ filter. Bahcall et al. (1997) noted that for diffraction-limited HST observations through a broadband filter, the color of an observed star (effectively, the SED-weighted average of the monochromatic PSFs) can significantly affect the quality of quasar point source subtractions if the quasar and star SEDs are poorly matched. To first order, redder objects will have a slightly broader PSF (e.g., as measured by FWHM), and bluer objects slightly narrower. Since we do not have detailed information about the shape of the quasar or star SEDs through the filter, and because of the degeneracy with focus, we simply leave the choice of PSF as a free parameter during the fitting procedure, which also allows matching of higher-order features (e.g., differences in diffraction spike patterns).

In the presence of detectable host galaxy flux, fitting only a single point source when attempting to remove the central quasar light tends to over-subtract the point source, especially for bright or centrally concentrated host galaxies. We, therefore, adopted a simultaneous fitting technique that models the point source and the underlying host galaxy flux distribution, approximating the latter with a Sérsic profile (following Mechtley et al. 2012). We stress that the actual morphologies of the host galaxies may differ from symmetric ellipsoids; the purpose of the Sérsic model is simply to provide a first-order approximation of the surface brightness distribution with a flexible parameterization.

We first attempted the two-component fit using the software GALFIT (Peng et al. 2002, 2010), currently the most widely used 2D surface brightness modeling software. However, this software employs several design decisions that make it less desirable for this particular problem (e.g., as noted in Section 6.2 of Peng et al. 2010).

First, GALFIT uses a standard Levenberg–Marquardt gradient descent method to perform least-squares minimization when fitting models. This method involves calculating a gradient image during each iteration to determine the parameter values to use for the subsequent iteration. Extremely compact models—e.g., a Sérsic profile with small effective radius (r_e) and large index (n)—have most of their gradient information contained within a single pixel, and so the Sérsic degrees of freedom are prone to fitting aberrant pixels (e.g., from PSF mismatch). This essentially creates a false minimum in parameter space, from which the gradient descent cannot escape.

A second, related problem is that GALFIT assumes that the supplied PSF is without error. Even without systematic PSF uncertainties—i.e., a PSF exactly matching the telescope focal history, spectral energy distribution of the quasar point source, etc.—the photon noise can be large enough to cause problems. While our PSF stars were selected to have higher S/N than the quasars, they exceed the quasar S/N by less than a factor of 10. This means that when performing the PSF subtraction, the PSF can contribute up to ≃30% of the per-pixel rms error.

Finally, the model is expected to have covariant parameters, such as the relative flux normalizations of the point source component and Sérsic component. These covariances are not quantified by GALFIT.

To address these problems, we developed our own Markov Chain Monte Carlo (MCMC) simultaneous fitting software, PSFMC (Mechtley 2014).¹⁰ As an MCMC parameter estimator, it addresses the first and third problems intrinsically, by allowing the user to provide prior probability distributions for the model parameters, and providing as an output product the full posterior probability distribution of model parameters given the observed data. The second problem is addressed by propagating a supplied PSF variance map during the convolution process.

Our software is built upon the PYMC Python module for Bayesian stochastic modeling (Patil et al. 2010). The computational book-keeping tasks—such as implementing the MCMC sampler, setting proposal distributions (see below), evaluating prior probabilities for a given set of parameters, and saving sample traces to disk—are handled by PYMC. We provide a framework that allows the user to simultaneously model an arbitrary number of model components (at this time, point sources, Sérsic profiles, and sky background). The free parameters for each component (e.g., position, total magnitude, Sérsic index, etc.) can either be supplied as a fixed numeric value or as an arbitrary prior probability distribution. An arbitrary number of samples can be drawn from the posterior distribution.

Samples are drawn from the posterior distribution using the Metropolis–Hastings algorithm (Metropolis et al. 1953; Hastings 1970), which accepts a proposed sample probabilistically based on the ratio of its posterior probability to that of the previous sample in the chain. Each time a proposed sample is drawn from the parameter space, PSFMC generates a model image of the intrinsic surface brightness distribution described by the parameters (hereafter, "raw model"). This raw model is then used to generate two further images—one convolved with the PSF ("convolved model"), and a model variance map.

The model variance map is simply the square of the raw model image convolved with the PSF variance map. The intensity, I_C(p), of a pixel p in the convolved model is given in Equation (1) below, where the summation is over all pixels q. I_R(q) is the raw model intensity, and $K(p,q)$ is the PSF weight for pixel p with the kernel centered at pixel q. If the kernel weights have associated variances ${\sigma }_{K}^{2}(p,q)$, then the variance propagates in the usual way, and the convolved variance of pixel p is given by Equation (2) below. This assumes that the off-diagonal covariance terms are zero (the noise between pixels is not correlated), which is almost true for the drizzle parameters used, and is a standard assumption in 2D surface brightness modeling

$$\begin{eqnarray}&&{I}_{C}(p)=\displaystyle \sum _{q}{I}_{R}(q)K(p,q),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\sigma }_{C}^{2}{(p)=\displaystyle \sum _{q}{I}_{R}(q)}^{2}{\sigma }_{K}^{2}(p,q).\end{eqnarray} \tag{ 2 }$$

The conditional probability of each observed pixel value given the model is then calculated, using a normal distribution with mean equal to the PSF-convolved model, and a variance equal to the summed convolved model and observed variances. The likelihood function is then the joint probability of these individual pixel probabilities, which is then multiplied by the prior probabilities of the parameter values to generate the sample's posterior probability. The sample is then accepted or rejected based on the Metropolis criterion.

Each quasar in our sample was modeled by two simultaneously fit components—a point source and a Sérsic profile. Although the software enables more complex multi-component fitting, the number of required samples (and thus computation time) increases steeply with additional parameters. Since we needed to model a total of 179 images (19 quasar images and 160 comparison galaxy images, see Section 4.1), this argued against iterative hand-crafting of more complex models for hosts that show multiple nuclei or other more complicated structures. Instead, we masked out other galaxies in each image with flux peaks that were separate from the central source (effectively, any galaxy with $\gtrsim 1\buildrel{\prime\prime}\over{.} 5$ separation), excluding them from the fit. This ensures that the Sérsic profile free parameters are used to fit the light distribution of the host galaxy (or galaxies) surrounding the quasar, rather than neighboring or foreground galaxies. This 2-component model is a computational compromise, with one Sérsic profile being more appropriate than none to avoid over-subtracting the point source. So we caution against over-interpreting the model parameter estimates. For instance, positional offsets between the point sources and the Sérsic profile centers should be interpreted as evidence for asymmetric flux distributions or multiple components in the host, rather than physical separation between the black hole and the center of its host.

The prior distributions adopted for the model parameters are summarized in Table 2. At z ≃ 2, the drizzled 0060 linear pixel scale corresponds to a physical size of 0.52 kpc. The ranges on the parameter priors were selected to model the entire range of values that are both physically reasonable and detectable. In particular, the Weibull distribution is used to approximate the observed distribution of Sérsic indexes (Yoon et al. 2011; Ryan et al. 2012).¹¹

Table 2.  Adopted Prior Distributions of Model Parameters

Model Parameter	Distribution	Value Range
Whole-image Parameters
Stellar PSF Image	Discrete Uniform	1–8 (list index)
Sky Background	Normal	$0\pm 0.01$ e−s−1
Point Source Component
X, Y Position	Normal	Centroid ± 4 pix
Total Magnitude	Uniform	$m{{}_{{\rm{H}}}}_{+1.5}^{-0.2}$ mag
Sérsic Component
X, Y Position	Normal	Centroid $\pm 8$ pix
Total Magnitude	Uniform	${m}_{{\rm{H}}}-26$ mag
Eff. Radius (Major)	Uniform	$1.0-15.0$ pix
Eff. Radius (Minor)	Uniform	$1.0-15.0$ pix
Index n	Weibull	$\alpha =1.5,\beta =4$
Position Angle	Circular Uniform	0°–180°

Note. Ranges are for intrinsic quantities, before convolution with the PSF. "Centroid" and m_H refer to the flux centroid and total $F160W$ magnitude of the (point-source dominated) quasar+host galaxy in the WFC3 image. Images are sky-subtracted during drizzling, but with some uncertainty, so the sky value is left as a free parameter.

Download table as:  ASCIITypeset image

Four chains were run for each model, each 100,000 samples with the first 50,000 discarded as a burn-in period, allowing the chains an opportunity to converge in parameter space before they are retained for analysis. Each chain is fully independent, with no parallel tempering of proposal distributions used. We use the Gelman–Rubin Potential Scale Reduction Factor to assess convergence ($\widehat{R}$, Gelman & Rubin 1992). This summary statistic compares the variance within individual chains to the variance among the chains to estimate how much sharper the posterior distribution for a parameter could be made by running for more samples. To consider a parameter's posterior estimate converged, we require $\widehat{R}\lt 1.05$, i.e., the estimated potential reduction in the scale of a parameter's univariate marginalized posterior distribution is less than 5%. Samplers whose chains had not converged (one or more parameters had $\widehat{R}\geqslant 1.05$) after 100,000 samples were continued for another 50,000 samples to increase the burn-in time, up to two additional times. Objects requiring longer convergence times are generally those where the host galaxies are only marginally detected. In all cases the point source parameters (position and magnitude) were well determined and met the above convergence criterion. In four of the 19 models, one or more Sérsic parameters still had $\widehat{R}\geqslant 1.05$. Since the goal of the modeling is point source subtraction, and single Sérsic models are a contrived simplification of the true flux distributions, we consider these converged for the purpose of our experiment.

The result of the MCMC fitting process is a collection of samples representing the posterior probability distribution of the free parameters in Table 2, given the observed data. The PSFMC software uses this to create model images—both before and after convolution with the PSF—and an image with all point source components subtracted. These images can either be made from the single maximum a posteriori (MAP) sample, or posterior-weighted. Posterior-weighted here means that an average image is made from the model images of all the retained MCMC samples (hence, weighted by the posterior), rather than the single image from the MAP sample. The purpose of this is primarily diagnostic—for models with well constrained parameter estimates, the posterior-weighted and MAP images are almost identical. For models with poorly constrained parameters or multiple posterior modes (e.g., SDSS J135851.73+540805.3), these will be visible in the posterior-weighted images, but not the MAP images, since MAP represent parameter estimates from only a single sample, rather than the full distribution of parameter values. Figure 1 shows the observed quasar, the posterior-weighted model image before PSF convolution, and the image of the host galaxy after subtracting only the point source component of the model from the observed data, both with original sampling and smoothed by a 2 × 2 pixel Gaussian.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Several cosmetic features in these images are worth discussing. In some of the raw model images ("Posterior Model" panels), a central dot is visible. This is the location of the point source before convolution with the PSF. The point source-subtracted images for several quasars (e.g., SDSS J113820.35+565652.8 and SDSS J124949.65+593216.9) show negative cores with a ring-like residual structure at $\simeq 0\buildrel{\prime\prime}\over{.} 22$ radius (the location of the first maximum in the Airy disk). Others (e.g., SDSS J102719.13+584114.3 and SDSS J155447.85+194502.7) show strong positive core residuals (most apparent in the smoothed images). These features are related to PSF focus and/or SED mismatch rather than significant point source over- or under-subtraction, and also appear when subtracting stars from other stars. The per-pixel S/N is lowest in the center due to the point source shot noise, and these mismatch structures are faint compared to the total point source signal, so they do not significantly affect the fitting and have per-pixel S/N $\lt \,1$ in the point source-subtracted images.

Having samples from the full posterior distribution also allows us to examine parameter covariance, as discussed in Section 3.2. The strongest of these covariances is between the total magnitude of the point source component and the total magnitude of the Sérsic component (though we note that the point source magnitude is still determined to within 0.01 mag on average, and the Sérsic magnitude to 0.1 mag). This is expected, since their combined flux must in some sense add up to the observed flux in the WFC3 image, but since some of the total flux is buried in the sky noise the exact form of the covariance function varies from object to object (arguing against removing one of the magnitude free parameters from the model). Figure 2 plots the posterior distribution for a typical object, marginalizing over each pair of parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For non-morphological analysis (see Section 6.5 below), we measure revised $F160W$ magnitudes for all objects, rather than using the Sérsic approximations resulting from the MCMC process. This is primarily important for deriving upper limits for the non-detections, since in those cases the Sérsic approximation often fits PSF mismatch structures rather than host flux (see e.g., the posterior model image for SDSS J135851.73+540805.3). Fluxes for the 16 detected galaxies, on the other hand, are identical to their Sérsic fits, within the 1-sigma magnitude errors. We perform isophotal photometry on the point source-subtracted images, using a threshold of S/N $=\,0.9$ times the rms background noise to determine the lowest isophote. In the galaxy centers, the high shot noise—including contributions from the removed quasar point source and the model PSF—results in few pixels having significant flux values (whether positive or negative, i.e., S/N is lowest in the center), and are primarily due to PSF mismatch (see discussion in Section 3.4 above). We, therefore, use values from the model for pixels within a $0\buildrel{\prime\prime}\over{.} 6$ diameter circle surrounding the point source location. For galaxies with total S/N $\lt \,2.0$ (the three non-detections), we report the $2\sigma$ upper limit. The difference compared to the Sérsic fit magnitudes is ${\rm{\Delta }}m\lt 0.5$ mag for the detected objects, but for non-detections the upper limits are $1.0\mbox{--}1.5$ mag brighter than the values indicated by the fit. The observed $F160W$ magnitudes and corresponding V-band absolute magnitudes are reported in Table 3 in Appendix A. The host-to-nuclear luminosity ratio ${L}_{\mathrm{Host}}/{L}_{\mathrm{Nuc}}$ of our sample spans the range from $\lt 0.01$–0.16, with median 0.058 (${m}_{\mathrm{Nuc}}\mbox{--}{m}_{\mathrm{Host}}=-3.2$ mag). This is roughly the range predicted for quasars with $L/{L}_{\mathrm{Edd}}\gtrsim 0.1$ from SDSS spectral decompositions (which are sensitive only to host fractions $\gt 0.1$, Vanden Berk et al. 2006).

Table 3.  Quasar Properties in Consensus Rank Order

Quasar	mNuc	mHost	MV	MBH	$L/{L}_{\mathrm{Edd}}$	Rank	SEM	ND Count
(SDSS J)	(mag)	(mag)	(mag)	$\mathrm{log}({M}_{\odot })$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
124949.65+593216.9	18.22	21.48 ± 0.27	−23.71	9.5	0.35	6	3.8	0
131535.42+253643.9	18.64	20.76 ± 0.10	−24.18	9.1	0.71	12	7.1	0
155447.85+194502.7	18.94	21.92 ± 0.16	−23.35	9.7	0.14	17	8.6	0
215954.45−002150.1	16.73	20.58 ± 0.36	−24.44	9.8	0.74	19	6.9	0
232300.06+151002.4	18.19	21.24 ± 0.21	−23.82	9.7	0.19	21	7.1	0
215006.72+120620.6	18.72	21.13 ± 0.11	−23.94	9.5	0.29	24	11.1	0
145645.53+110142.6	18.21	20.94 ± 0.14	−24.19	9.2	0.65	28	7.9	0
120305.42+481313.1	18.24	20.98 ± 0.16	−24.09	9.4	0.27	33	8.8	0
143645.80+633637.9	17.00	20.41 ± 0.29	−24.81	10.1	0.22	37	12.7	1
220811.62−083235.1	18.53	21.92 ± 0.23	−23.01	9.2	0.63	38	8.7	0
131501.14+533314.1	18.51	22.26 ± 0.37	−22.67	10.0	0.10	41	9.2	0
094737.70+110843.3	19.07	21.39 ± 0.11	−23.51	9.7	0.11	50	14.3	0
081518.99+103711.5	18.97	21.20 ± 0.13	−23.92	9.3	0.49	66	8.0	0
082510.09+031801.4	17.74	(21.88)	(−23.28)	9.7	0.28	69	10.6	1
113820.35+565652.8	18.14	21.95 ± 0.41	−22.97	9.5	0.30	71	11.1	0
085117.41+301838.7	18.88	22.74 ± 0.34	−22.18	9.0	0.65	72	9.6	0
102719.13+584114.3	18.68	(22.8)	(−22.33)	9.3	0.57	85	11.6	4
123011.84+401442.9	17.44	19.56 ± 0.11	−25.63	9.9	0.28	88	11.2	2
135851.73+540805.3	18.25	(22.25)	(−22.97)	9.4	0.54	95	5.5	9

Note. Quasar properties, in consensus rank order. Magnitude values in parentheses indicate $2\sigma$ upper limits. Column 1: quasar name, giving the full sexagesimal coordinates; Column 2: observed magnitude of the subtracted nuclear point source; Column 3: observed magnitude of the host galaxy; Column 4: rest-frame V-band absolute magnitude of the host galaxy; Column 5: black hole mass (see Table 1); Column 6: Eddington specific accretion ratio (see Table 1); Column 7: consensus rank; Column 8: standard error on the mean of consensus rank; Column 9: non-detection count, i.e., number of individuals who flagged the host as a non-detection.

Download table as:  ASCIITypeset image

Table 4.  Inactive Galaxy Properties in Consensus Rank Order

Field	ID	R.A.	Decl.	Redshift	${m}_{F160W}$	M*	AV	Rank	SEM	Use Count	ND Count
		(degree)	(degree)	z	(mag)	$\mathrm{log}({M}_{\odot })$	(mag)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
GOODS-N	16827	189.297348	+62.225300	1.97	21.47	11.21	1.3	0	1.0	1	0
AEGIS	38187	214.75116	+52.830135	2.13	21.52	11.51	2.5	1	15.0	1	1
GOODS-N	3253	189.007736	+62.150768	1.98	22.00	10.36	0.4	2	16.5	1	1
COSMOS	28406	150.191299	+2.480364	2.02	21.97	10.49	0.7	3	2.3	1	0
UDS	35887	34.379311	−5.156958	1.93	21.32	10.95	0.9	4	2.3	1	0
COSMOS	11363	150.119614	+2.295948	2.12	21.31	11.07	1.1	5	2.5	3	2
AEGIS	13874	215.097092	+52.977772	1.98	21.69	10.27	0.3	7	6.0	1	0
AEGIS	38157	214.643967	+52.751797	1.83	21.80	10.97	2.2	8	3.2	1	0
AEGIS	41042	214.59761	+52.739277	1.82	21.72	10.76	0.5	9	3.3	1	0
AEGIS	3311	215.045425	+52.894630	1.85	21.15	11.33	0.6	10	3.5	2	0
AEGIS	6310	214.780258	+52.719028	1.89	21.32	11.12	0.6	11	7.0	2	1
COSMOS	10989	150.077515	+2.292315	1.85	21.85	11.00	0.5	13	13.1	1	0
AEGIS	11930	215.078857	+52.956783	1.75	21.99	9.90	005	14	6.9	1	0
COSMOS	28565	150.139343	+2.481772	2.07	21.88	11.4	0.6	15	19.5	1	2
AEGIS	38652	215.093567	+53.074688	2.22	22.50	10.43	0.7	16	22.4	1	4
AEGIS	22887	214.666656	+52.707596	1.89	22.49	11.11	1.3	18	9.8	1	0
GOODS-S	30274	53.131084	−27.773108	2.24	21.31	11.16	1.2	20	5.0	3	0
UDS	25091	34.308601	−5.192824	1.80	21.34	10.95	0.3	22	5.7	1	0
GOODS-N	30283	189.444336	+62.294060	2.10	20.61	11.32	1.1	23	5.9	5	17
COSMOS	10592	150.079071	+2.288243	1.76	21.73	11.02	0.3	25	9.1	1	0
GOODS-N	10657	189.256592	+62.196178	2.22	21.58	11.68	1.9	26	5.5	2	0
AEGIS	15871	215.030258	+52.937485	2.17	21.91	11.07	1.2	27	5.2	2	0
GOODS-N	20709	189.464142	+62.244041	1.80	21.64	11.23	0.7	29	5.5	1	0
COSMOS	4536	150.13826	+2.225018	1.80	21.89	11.17	2.4	30	6.3	1	0
AEGIS	36755	215.086044	+53.060520	1.94	21.63	10.52	0.6	31	10.8	1	0
COSMOS	2716	150.177567	+2.208040	1.94	21.79	10.81	1.0	32	10.2	1	0
AEGIS	29591	214.799545	+52.829071	1.77	21.96	10.45	1.1	34	8.4	2	0
GOODS-S	4568	53.164497	−27.890385	2.13	21.37	10.66	0.3	35	9.8	1	0
AEGIS	5745	215.104919	+52.949425	1.76	22.48	10.63	0.6	36	16.5	1	2
UDS	922	34.360863	−5.275831	1.79	21.93	11.24	0.7	39	7.1	1	0
COSMOS	17406	150.158478	+2.357966	1.79	21.34	11.10	1.2	40	8.1	2	0
COSMOS	13083	150.0961	+2.313470	2.01	21.76	11.05	0.6	42	5.8	2	1
COSMOS	7884	150.065628	+2.260996	2.05	21.27	11.44	0.6	43	8.4	3	0
GOODS-S	37243	53.098831	−27.736591	1.76	21.42	10.79	1.2	44	12.8	1	0
GOODS-S	36095	53.162609	−27.739000	1.96	21.18	11.45	0.6	45	9.2	3	0
AEGIS	16065	214.979965	+52.902481	2.19	22.51	11.14	2.6	46	9.5	1	0
AEGIS	24059	215.017929	+52.962994	1.86	22.54	11.07	1.0	47	10.0	1	0
AEGIS	25526	215.101151	+53.026985	1.77	21.94	10.75	0.6	48	12.0	3	12
GOODS-S	10562	53.266895	−27.861710	1.90	21.76	10.80	0.7	49	8.0	1	0
UDS	14996	34.278858	−5.226638	2.07	21.74	11.20	1.3	51	16.3	1	1
AEGIS	29790	215.081818	+53.030506	1.92	22.46	9.74	0.2	52	15.5	1	6
AEGIS	24333	215.137421	+53.046440	1.97	21.56	11.47	0.9	53	9.9	2	0
AEGIS	25734	214.837219	+52.840817	1.88	22.31	10.36	0.9	54	16.6	1	2
GOODS-N	2694	189.208069	+62.146435	2.10	21.95	10.84	0.4	55	14.1	2	3
UDS	36685	34.289249	−5.153751	1.95	21.80	11.29	1.0	56	10.1	1	0
UDS	30475	34.544548	−5.174858	1.75	21.24	10.81	0.7	57	9.3	4	9
AEGIS	38153	214.707108	+52.797298	1.90	21.83	10.91	2.4	58	14.5	2	0
AEGIS	23768	214.780533	+52.792133	1.88	21.84	10.82	1.1	59	9.4	1	0
AEGIS	20962	215.126419	+53.027435	2.18	22.51	10.41	1.0	60	9.5	1	0
AEGIS	18678	214.881332	+52.843845	1.99	22.34	10.49	1.0	61	13.8	1	2
AEGIS	11730	215.184616	+53.031445	1.98	22.48	11.1	0.6	62	10.3	1	1
COSMOS	11136	150.109085	+2.294313	1.88	21.89	10.90	2.0	63	10.1	2	0
COSMOS	21413	150.088943	+2.399108	1.78	21.77	10.97	1.2	64	7.5	2	0
UDS	23692	34.363182	−5.199362	2.16	21.04	11.60	1.1	65	8.0	5	3
COSMOS	7411	150.177017	+2.255223	1.79	21.50	11.00	0.6	67	11.3	1	0
UDS	4721	34.2393	−5.263094	1.92	21.38	11.4	0.0	68	10.0	1	0
AEGIS	296	215.247437	+53.020676	1.98	21.98	10.93	0.9	70	12.3	1	0
UDS	25630	34.372494	−5.191559	1.78	20.79	11.26	1.3	73	10.5	4	5
AEGIS	2918	215.168991	+52.980843	2.07	21.32	11.53	1.0	74	9.5	3	0
AEGIS	17644	214.906448	+52.857258	1.79	21.81	10.99	1.2	75	10.1	1	1
COSMOS	25581	150.154617	+2.444321	1.86	21.05	11.31	0.2	76	9.1	5	9
AEGIS	36574	215.094772	+53.066292	2.07	21.62	11.19	0.1	77	8.5	2	2
COSMOS	19071	150.087845	+2.373441	1.79	21.76	10.82	0.3	78	8.5	1	0
COSMOS	11494	150.073929	+2.297983	2.05	20.84	11.57	0.9	79	5.3	4	14
AEGIS	2016	215.226395	+53.017464	1.98	22.52	11.00	0.7	80	9.1	1	0
AEGIS	26649	215.042541	+52.989807	1.87	21.71	11.13	1.3	81	6.4	2	0
COSMOS	16419	150.095642	+2.350081	2.09	21.09	11.29	0.0	82	8.5	5	22
UDS	4701	34.273331	−5.262073	2.11	21.11	11.18	0.4	83	8.3	2	3
COSMOS	20983	150.108505	+2.393802	1.97	21.66	11.07	0.9	84	9.5	1	0
COSMOS	8443	150.079056	+2.266573	1.82	21.11	10.95	0.7	86	6.9	7	19
GOODS-N	2295	189.11438	+62.142494	2.04	21.78	11.03	0.5	87	6.9	1	8
UDS	32892	34.389614	−5.168086	1.90	21.24	10.98	0.3	89	6.9	2	4
AEGIS	9128	214.976212	+52.872486	2.13	21.87	11.32	0.5	90	9.5	1	3
GOODS-N	33780	189.202026	+62.317154	1.86	21.41	11.29	1.1	91	7.8	2	0
UDS	20529	34.549076	−5.208462	1.77	21.00	11.26	0.7	92	5.7	5	9
GOODS-N	23187	189.514648	+62.255245	1.77	21.1	11.08	0.5	93	7.6	6	9
UDS	22802	34.446915	−5.200699	1.77	21.05	10.99	0.0	94	8.6	5	9
GOODS-N	17735	189.060898	+62.228977	1.84	21.48	10.94	0.4	96	8.0	1	1
UDS	19405	34.544777	−5.211998	1.86	21.75	10.92	0.4	97	8.3	2	9
AEGIS	12417	215.052643	+52.940708	1.82	22.46	10.94	0.3	98	4.5	1	8
AEGIS	35002	215.093201	+53.058929	1.79	21.92	11.07	2.2	99	6.8	3	2
AEGIS	12288	215.131424	+52.995796	1.98	21.86	10.83	0.4	100	9.9	1	3
GOODS-S	42501	53.197029	−27.712065	1.82	21.70	11.00	0.6	101	6.4	1	7
COSMOS	17255	150.133347	+2.355615	1.76	21.71	10.82	0.3	102	8.8	1	3

Note. Inactive galaxy properties, in consensus rank order. Column 1: CANDELS field name; Column 2: 3D-HST catalog ID number; Column 3: R.A.; Column 4: decl.; Column 5: Redshift (3D-HST z_peak); Column 6: observed magnitude in the $F160W$ filter; Column 7: stellar mass (from 3D-HST SED fit); Column 8: rest-frame V-band internal dust attenuation (from 3D-HST SED fit); Column 9: consensus rank; Column 10: standard error on the mean of consensus rank; Column 11: number of quasars for which the galaxy was selected as comparison; Column 12: non-detection count, i.e., number of individuals who flagged the host as a non-detection (noting the maximum is Use Count × ten experts).

Download table as:  ASCIITypeset images: 1 2

As discussed in Section 1, any study that wishes to make a definitive statement about particular AGN triggering mechanisms needs to compare the AGN hosts to a matched sample of inactive galaxies. This matching, in principle, should be done for all properties except accretion rate (AGN versus non-accreting) and the property we wish to test for (morphology). Our quasars were selected by SMBH mass, but this is unfortunately not a directly measurable quantity for large samples of inactive z = 2 galaxies in HST extragalactic survey fields. Stellar masses (M_*) are, on the other hand, fairly well constrained in such fields from multi-wavelength SED fitting, when the available data cover the entire rest-frame ultraviolet to rest-frame near-infrared range (with the standard caveats of stellar population synthesis modeling, including an assumed universal initial mass function). The question is then how to select galaxies whose stellar masses match the distribution of stellar masses for the quasar hosts. One approach is to assume a certain ${M}_{\mathrm{BH}}\mbox{--}{M}_{* }$ relation (e.g., Häring & Rix 2004; Kormendy & Ho 2013), taking into account scatter in the relation, possible redshift evolution, and biases introduced by our selection function, and then select galaxies in the resulting stellar mass range. The second approach is to match samples using some proxy for M_* that is directly observable in both data sets, namely the total $F160W$ magnitude. We opted for this second option, since it requires no assumptions regarding SMBH scaling relations and still captures the full range of galaxies that might reasonably host the high-mass SMBHs. We discuss scaling relations in more detail in Section 6.5.

We obtained $F160W$ images for our comparison sample from the CANDELS HST multi-cycle treasury program (Grogin et al. 2011; Koekemoer et al. 2011), selected using the redshift and photometry catalogs from the 3D-HST program (v4.1, Brammer et al. 2012; Skelton et al. 2014). The five CANDELS treasury fields represent the best comparison data set since they image a moderately wide area in $F160W$ (0.22 deg²), and contain a wealth of ancillary data at other wavelengths to ensure accurate redshift determinations, stellar masses, and AGN identifications. We first selected all galaxies in the redshift range z = 1.8–2.2. Because the CANDELS fields also have deep $3.6\mbox{--}8.0\,\mu {\rm{m}}$ imaging from the Spitzer Space Telescope, which samples the red side of the Balmer/4000 Å break at z = 2, the 3D-HST photometric redshifts are sufficiently accurate (${\rm{\Delta }}z/(1+z)\lt 0.1$ at z = 2, Skelton et al. 2014). We next exclude probable AGN from the sample by cross-matching the z ≃ 2 galaxies with X-ray catalogs (Subaru/XMM-Newton Deep Survey for UDS, Ueda et al. 2008; Chandra Source Catalog for the other fields, Evans et al. 2010). The number of X-ray sources removed from the z ≃ 2 sample in each field were: AEGIS: 9, COSMOS: 5, GOODS-N: 4, GOODS-S: 6, UDS: 1. We examined each remaining object and excluded any spurious detections (stars and diffraction spikes). This left us with a parent sample of 1123 z ≃ 2 galaxies with ${m}_{F160W}\lt 23$ mag, and 150 with ${m}_{F160W}\lt 22$ mag.

We then drew ten samples from the posterior distribution of Sérsic component magnitude for each of the 16 detected quasar hosts¹² , and selected the z ≃ 2 galaxies which most closely matched these samples in $F160W$ magnitude. The quasar hosts are luminous (median ${m}_{F160W}=21.36$ mag), so there are not enough objects in CANDELS to ensure each host has ten unique magnitude-matched comparison galaxies. We, therefore, allowed galaxies to be re-used as comparisons for more than one quasar (i.e., drawn galaxies were replaced before selecting galaxies matched to the next quasar). In total, 48 comparison galaxies were used once, 18 were used twice, and 18 three or more times (with a maximum of seven times). This gave us ten redshift- and magnitude-matched inactive comparison galaxy images for each quasar, or a total of 160 images of 84 unique galaxies. We were originally prepared to either down-select to fewer comparison galaxies per quasar or employ weighting in analyses, but decided this was unnecessary since the weighted and unweighted stellar mass distributions were not significantly different (see Section 6.5 below). Our analyses simply consider the sample of 84 comparison galaxies without re-use, unless otherwise stated (though see Appendix A for some consistency checks that were enabled by this re-use).

To allow a fair comparison between quasar hosts and inactive galaxies, we needed to ensure that merger signatures were equally detectable in both sets of images. The CANDELS $F160W$ images are deeper than those in our SNAP program, and the quasar images retain systematic residual patterns from slight PSF mismatch, despite our best efforts at PSF matching (see Section 3.1). We, therefore, added synthetic point sources to the comparison galaxy images and performed the same MCMC point source subtraction procedure used for quasars.

For galaxies that were used as comparisons for multiple quasars (Section 4.1 above), the images were first transformed with a unique sequence of mirror flips and 90° rotations, so that no two images of a galaxy were exactly the same. For each comparison galaxy image, we then randomly selected one of the eight PSF stars (see Section 3.1) to act as the synthetic point source. This was then scaled to match the magnitude of the subtracted point source for the galaxy's corresponding quasar—recalling that ten comparison galaxies were selected to match each quasar host—and inserted into the image centered at the galaxy flux centroid. We then measured the sky noise in the resulting image (the point source having contributed some noise), and added noise if necessary to match the sky noise of the quasar image. This was approximated as an uncorrelated Gaussian field rather than correlated Poisson noise like the real quasar images, but this added field is not the dominant source of noise in the final images. Of the 160 images, 116 required no added noise, with the other 44 requiring at most $\sigma =0.0047$ e⁻s⁻¹ additional noise, corresponding to 29% of the final variance after point source subtraction. Weight images were produced by summing the individual variance sources for these images: the CANDELS (sky) variance map, the CANDELS object signal (approximated as the image in electron units), the scaled point source variance map, and the variance of the added Gaussian field. We then performed point source subtraction using the same MCMC parameter estimation procedure described in Section 3.3, supplying only the remaining seven stars as PSF models, i.e., the star used for each galaxy's synthetic point source was excluded. The resulting galaxy images then contained the same structural residuals from point source subtraction as the quasar host images.

After preparing the comparison sample, we now had images for a joint sample of 179 quasar hosts and inactive galaxies with similar noise properties and residuals from point source removal. Rather than classifying the galaxies into bins of merger strength or interaction stage (e.g., Cisternas et al. 2011; Kocevski et al. 2012), we instead had ten galaxy morphology experts¹³ rank the 179 galaxies by morphological evidence of distortion due to strong gravitational interactions. Since ranking is a relative measure of interaction strength, it gives a natural way to avoid personal statistical biases (e.g., what an individual considers a major versus a minor merger), and allows us to examine how our conclusions change as a function of where exactly the final distinction between merger and non-merger is drawn. While ranking the galaxies, these experts were also asked to note any instances where the galaxy was undetected (i.e., was below the noise limit or completely obscured by image artifacts from the point source subtraction).

We then combined the ranked lists to establish a single consensus sequence of galaxies from the most to the least distorted. This is more difficult than might naively be assumed, since it is equivalent to a ranked voting system and results from social choice theory apply. In particular, Arrow's Impossibility Theorem (Arrow 1950) states that there is no way to combine such ranked lists such that the consensus sequence has three desirable properties: (1) non-dictatorship, such that each voter receives equal weight, (2) unanimity, such that if all voters agree option $X\gt Y$, then the consensus sequence also has $X\gt Y$, and (3) independence of irrelevant alternatives, such that the consensus preference between X and Y depends only upon the individual preferences between X and Y, and not relationships between other options. Several techniques are available for relaxing one or more of these properties—we chose to relax non-dictatorship by allowing de-weighting of individual votes, clipping those that significantly disagreed with the majority for each galaxy.

The individual lists were merged by calculating the mean rank for each galaxy and its associated variance. We then clipped any individual ranks that were more than $2\sigma$ from the mean rank for each galaxy. This excluded 55 of the 1790 individual rankings, roughly evenly distributed between rankers (between zero and 11 individual ranks clipped for each ranker). A total of 237 individual rankings were additionally flagged as non-detections. For comparison galaxies that were included multiple times (see Section 4), we used the inverse variance-weighted mean of the individual rankings as the final rank. If half or more of the individual rankings flagged a galaxy as non-detected, we excluded it from further analysis (as noted in the merger statistics below). Images of all the galaxies in consensus rank order are included in Appendix A.

As with all studies based on visual inspection, individual galaxies show more or less evidence for major mergers. We made this explicit by ranking them on a continuum from the most to least evidence. An alternative to merger fractions, then, is to consider the two samples' distributions within this ranked continuum. These distributions are formally indistinguishable under either a 2-sample Kolmogorov–Smirnov test ($D=0.21,p=0.45$) or 2-sample Anderson–Darling test (${A}^{2}=-0.38,p=0.52$ ). The notion of categorizing galaxies as mergers or non-mergers is convenient for some of the discussion in Section 6, so we discuss the derivation of a merger fraction for each sample below.

To estimate the quasar host and inactive galaxy merger fractions, we visually inspected the consensus sequence and selected a particular cutoff rank, below which we could no longer find any clear merger signatures. We selected rank 30 for this cutoff (corresponding to 32% of the final sample of all quasars and inactive galaxies, minus the non-detections), but note that the particular choice of rank may differ among individual experts. In Appendix A we explore in detail how this choice affects the inferred merger fractions. We only note here that for any reasonable choice of cutoff rank, the qualitative aspects and conclusions below do not change.

With a merger/non-merger cutoff rank of 30, we identified seven quasar hosts with evidence of major mergers, 11 with no evidence of major mergers, and one indeterminate (more than half of individuals flagged them as non-detections). For the inactive galaxies, 24 had evidence of major mergers, 56 had no evidence of major mergers, and four were indeterminate. The probability distribution describing the inferred merger fraction given the finite number of experimental trials is the beta distribution. The PDF for the beta distribution is given by Equation (3):

$$\begin{eqnarray}&&P{(x;a,b)=\displaystyle \frac{(a+b+1)!}{a!b!}{x}^{a}(1-x)}^{b}.\end{eqnarray} \tag{ 3 }$$

Here a and b are integers denoting the number of mergers and non-mergers, respectively. The resulting probability distributions for both quasar hosts and inactive galaxies are shown in Figure 3. The modes or peaks of the distributions are simply the usual merger fractions, $a/(a+b)$. The inferred parent population merger fractions, with 68% confidence intervals, are thus ${f}_{{\rm{m}},\mathrm{qso}}=0.39\pm 0.11$ for quasar hosts, and ${f}_{{\rm{m}},\mathrm{gal}}=0.30\ \pm 0.05$ for inactive galaxies.

Zoom In Zoom Out Reset image size

Although we use the term "merger fraction," we emphasize that the distorted fractions above should not be used as absolute merger fractions to compare to other studies of high-mass galaxies. We, as with other studies, have invariably missed a few real mergers and possibly mis-identified some non-mergers due to noise and residual patterns associated with the point source subtractions (see discussion in Appendix B). Rather, they should only be interpreted relative to each other, since both samples have the same observational limitations. The probability distribution for the quasar merger fraction in Figure 3 rules out the extreme scenario immediately. The merger fraction is not consistent with values near 1 (99.7% confidence interval: $0.13\mbox{--}0.72$), as might be expected if every massive quasar were growing due to ongoing merger activity.

The quasar host PDF peaks at a slightly higher merger fraction than the inactive galaxy PDF, corresponding to a merger fraction enhancement of ${f}_{{\rm{m}},\mathrm{qso}}/{f}_{{\rm{m}},\mathrm{gal}}=1.3$. The enhancement is not significant, however. Given two random variables X and Y, the probability that $X\gt Y$ is the integral of the joint PDF over the region where this inequality holds. For the two distributions in Figure 3, this probability is $P({f}_{{\rm{m}},\mathrm{qso}}\gt {f}_{{\rm{m}},\mathrm{gal}})=0.78$, or $0.78\sigma$. We can invert this to examine the statistical sensitivity of our observations—given the fixed quasar sample size, and assuming the observed PDF for the inactive galaxies, we can calculate the minimum merger fraction that would have resulted in a significant ($\gt 3\sigma$) signal. For 18 quasars, this would have required 13 mergers, for a merger fraction of ${f}_{{\rm{m}},\mathrm{qso}}\geqslant 0.72$ or enhancement of ${f}_{{\rm{m}},\mathrm{qso}}/{f}_{{\rm{m}},\mathrm{gal}}\geqslant 2.4$. Correspondingly, if the true intrinsic distortion fractions are 0.30 and 0.39, we can ask how much larger the sample would need to be to detect this enhancement signal at a significant level. Both the quasar and comparison galaxy sample sizes would need to be increased by a factor of 7.6 to detect an enhancement signal at $2\sigma$, and a factor of 17.2 for $3\sigma$. This underscores the statistical insignificance of the observed enhancement, and the need for cautious interpretation when dealing with beta-distributed quantities inferred from small samples.

Although the observed enhancement is not significant, the data are still consistent with the quasars having either a slightly enhanced merger fraction, or no enhancement. The data are not, however, consistent with the quasars having a large merger fraction enhancement. If the observed enhancement is real (i.e., if it were still present with much larger samples), it could indicate that mergers are simply one of several possible AGN triggering mechanisms, rather than the dominant or only mechanism.

One possible caveat is that of significant time lags between the appearance of merger signatures and the onset of quasar activity. That is, if gas disturbed or provided by a merging companion takes long enough to reach the SMBH (i.e., several times the dynamical timescale), or is first reprocessed via an episode of star formation, morphological signatures of a merger may no longer be observable. Observational evidence for such time lags is necessarily indirect, and the most successful studies have examined growth regimes very different from those in our study—e.g., Wild et al. (2010), who studied AGN 2.5 orders of magnitude less massive than ours at z = 0.01–0.07, and inferred a time lag of 250 Myr between peak starburst activity and peak AGN accretion. Theoretical models of merger-induced AGN activity predict a delay of $\simeq 100\,{\rm{Myr}}$ between galaxy coalescence and the peak of quasar activity (Di Matteo et al. 2005; Hopkins et al. 2008). These estimates are at any rate shorter than the timescale over which morphological merger signatures are still visible (as high as $\gtrsim 1\,{\rm{Gyr}}$ for gas-rich mergers Lotz et al. 2010a, 2010b). For exceptionally long time lags, such that merger signatures are almost completely erased, the observational signatures would become essentially indistinguishable from VDIs (e.g., Bournaud et al. 2012; Trump et al. 2014). In such a case, the theoretical model might even be indistinguishable. That is, if gas transport times are several times longer than the dynamical time, is the merger responsible for the loss of angular momentum, or are secular processes like VDI simply acting upon the gas-rich merged galaxy? For the specific feeding mechanism we are testing—near-zero angular momentum gas provided directly by a major merger—the relevant timescale should be closer to the free-fall time and thus shorter than the timescale for merger signatures to be smoothed out.

We also examined the incidence of merger features as a function of black hole mass and accretion rate (Figure 4), with the caveat that our sample spans only about a factor of ten in each. No trend is observed in distortion rank as a function of black hole mass and thus no evidence that the higher-mass black holes are nearer the end of a current merger. Similarly, no trend is observed as a function of either luminosity or Eddington ratio, as might be expected if nuclear gas availability were systematically greater for merger-triggered quasars. We have planned a future study to address the question of merger triggering in high-Eddington systems specifically, as contrasted with low-Eddington AGN and inactive galaxies.

Zoom In Zoom Out Reset image size

Environmental diagnostics provide another test of merger or interaction triggering hypotheses. Small-scale environmental studies test a slightly different hypothesis compared to our study—essentially, that interactions trigger gas instabilities within a galaxy at early merger stages or in non-merging close encounters. Studies of small-scale clustering of quasars (Hennawi et al. 2006; Myers et al. 2008) found some evidence for an excess in the quasar spatial autocorrelation function at small (≃10 kpc) scales. They attributed this excess to gravitational interaction events triggering quasars in rich environments. Silverman et al. (2009) examined the AGN fraction in galaxies as a function of their local environmental density using a nearest-neighbors approach. They found that the hosts of these lower-luminosity AGN generally trace the same environments as star-forming galaxies (both processes requiring gas), with a preference toward under-dense regions for AGN hosts comparable in mass to our inactive sample (${M}_{* }\gtrsim {10}^{11}\,{M}_{\odot }$).

A related diagnostic for examining AGN triggering via gravitational interactions is the study of galaxies in kinematic pairs, i.e., those that are close in mass, spatial projection, and line-of-sight velocity, looking once again for an enhancement to the number of AGN in close pairs versus a field sample. Ellison et al. (2008) studied low-redshift ($z\lt 0.15$) kinematic pairs from SDSS (2402 galaxies in pairs, 69,583 in the field control sample), finding the (lower-luminosity) AGN fraction in pairs to be consistent with that of the control sample. Silverman et al. (2011) performed a similar study at z = 0.25–1.05 using the zCOSMOS spectroscopic sample (562 galaxies in pairs, 2726 control galaxies). They detected significant differences in the AGN fractions of three subsamples: close pairs (9.7%), wide pairs (6.7%), and isolated galaxies (3.8%). They estimated that 17.8% of moderate-luminosity nuclear activity is triggered during early-stage interactions, leaving the further ≃80% unaccounted for. Both the SDSS and zCOSMOS studies posited that later stages of major mergers may account for some of the missing triggers. Lackner et al. (2014) examined this by looking for galaxies with multiple nuclei at small separations (i.e., pairs no longer distinguishable as such in the low-resolution ground-based imaging), finding that the combination of wide pairs, close pairs, and late-stage mergers account for a total of 20% of AGN activity at z = 0.25–1.0. This is then consistent with the Cisternas et al. (2011) study, which found that morphologically identified major mergers are not a dominant trigger for low-to-moderate luminosity AGN. We have now further shown that such major mergers are not a dominant triggering mechanism for the high-mass quasars that dominate SMBH accretion at z = 2.

As mentioned in Section 1, studies of reddened and dust-obscured quasars have generally found very high merger fractions, in contrast with most studies of the unobscured AGN or quasar populations mentioned above. Various obscured or dust-rich quasar selection methods target slightly different spectral features. The first HST studies of red quasars targeted objects with far-infrared excesses (e.g., Canalizo & Stockton 2001), finding evidence for major mergers in eight of nine hosts. Although the authors argue that chance hosting is unlikely—i.e., the quasar and starburst activity are related—such objects make up a small fraction of low-redshift quasars, and so correspondingly represent a small contribution to the total triggering budget.

Zakamska et al. (2006) selected Type 2 quasars at z ≃ 0.2–0.4 from SDSS using the scattered emission-line luminosities as a proxy for total nuclear luminosity—i.e., objects expected to be similar to moderate-luminosity Type 1 quasars, but viewed from an angle where the circumnuclear dust torus obscures the direct AGN light, so that they are not distinct from an evolution standpoint. They found evidence of a major merger in only one of nine hosts, and evidence of tidal debris in three of nine, roughly in line with the lower merger fractions found in low-redshift Type 1 quasars (McLure et al. 1999; Dunlop 2001; Dunlop et al. 2003; Floyd et al. 2004).

More recently, considerable effort has gone into the study of dust-reddened Type 1 quasars selected from a combination of radio and near-IR data (F2M quasars, Glikman et al. 2007). These highly reddened quasars make up ≃10% of the luminous quasar population, are among the most intrinsically luminous at any given redshift (Glikman et al. 2012), and their host galaxies show a high incidence of mergers (merger fractions of ≃0.8–0.85, Urrutia et al. 2008; Glikman et al. 2015).

How do we reconcile these conflicting results with our current study and other unobscured quasar host studies? Reddening and obscuration can (at least) originate from non-evolutionary nuclear geometric effects (torus obscuration), evolutionary effects (a buried/blowout phase), and non-evolutionary host geometric effects (i.e., quasars that happen to be hosted in asymmetric, dust-rich, ULIRG-like galaxies). Any given red selection method may pick up some combination of the three. Type 2 selection like that of Zakamska et al. (2006) purports to select only for torus effects, and the similarity of those hosts to Type 1 hosts seems to support this. For red samples with high host merger fractions, additional evidence is needed to determine whether the mergers are an evolutionary feature or a host sub-population feature. The F2M quasars are the best bet for quasars triggered by early-stage galaxy mergers (see especially discussions in Glikman et al. 2012; Urrutia et al. 2012), but the universality of such triggering remains unclear. We argue in Section 6.2 above that, assuming the usual quasar lifetime estimates, we would see far more evidence for mergers in our sample if merger triggering were universal and all (or even most) quasar hosts begin as train-wreck mergers like the F2M hosts.

The distribution of the inactive galaxy stellar masses taken from the 3D-HST catalogs is plotted in Figure 5. Two histograms are shown: one where each galaxy is counted only once, and one where each is weighted by the number of times it was used for comparison to a quasar host. The median stellar mass of the unweighted comparison sample is $\mathrm{log}({M}_{* }/{M}_{\odot })=11.0$. The median stellar mass of the weighted sample (i.e., with the distribution of $F160W$ magnitudes matched to the quasar hosts) is $\mathrm{log}({M}_{* }/{M}_{\odot })=11.1$. We note that CANDELS/3D-HST are complete to much lower masses. There is no significant trend in stellar mass as a function of distortion rank, consistent with previous studies finding no significant mass-dependence of the major merger rate over similar stellar mass ranges (e.g., Xu et al. 2012).

Zoom In Zoom Out Reset image size

We can obtain crude estimates of the quasar host galaxy stellar masses from their $F160W$ magnitudes by adopting some assumptions about their stellar populations. Since we sample the SED at a longer wavelength than the Balmer/4000 Å break, the light comes mostly from older stars that account for the bulk of stellar mass. The general procedure is to calculate a galaxy's luminosity from some bandpass, then multiply by a stellar mass-to-light ratio (M/L). This assumes a particular SED, which may differ from object to object. Alternatively, we could use the empirical relation (with scatter) between $F160W$ magnitude and stellar mass from the Skelton et al. (2014) catalog, which implicitly encodes the full range of observed SEDs or M/L ratios for galaxies at a given redshift and stellar mass. The quasar hosts and inactive galaxies show a similar range of morphologies, so we take this approach.

Since this is the same population of galaxies from which the comparison sample was drawn, and since the quasar and inactive galaxy samples have the same luminosity distribution by construction, this should result in a mass distribution roughly the same as the inactive galaxies. Indeed, the quasar hosts have a stellar mass distribution of $\mathrm{log}({M}_{* }/{M}_{\odot })=11.2\pm 0.4$, but with an uncertainty on the mean value of 0.4 dex, due to the intrinsic scatter in the ${m}_{F160W}\mbox{--}{M}_{* }$ relation (in turn due to the physical range in SEDs for a fixed magnitude). There is no statistically significant difference between the distributions of stellar mass for the mergers and non-mergers, as diagnosed by either a 2-sample Kolmogorov–Smirnov test ($D=0.33,p=0.68$) or 2-sample Anderson–Darling test (${A}^{2}=0.046,p=0.33$ ). The average ${\text{}}M/L\simeq 0.5$ is consistent with relatively young stellar populations found in previous quasar host studies (e.g., Jahnke et al. 2004a, 2004b; Sánchez et al. 2004).

With stellar mass estimates for the quasar hosts, we can now compare them to the local ${M}_{\mathrm{BH}}\mbox{--}{M}_{* }$ relation. Figure 6 shows the z = 2 quasar hosts alongside the local relation derived by Kormendy & Ho (2013). The error bars include contributions from the scatter in both the virial black hole mass calibration (0.3 dex) and the $F160W$ magnitude-stellar mass relation (0.4 dex, encoding the range of SEDs for a galaxy with a given magnitude). The median black hole-to-stellar mass ratio is ${\rm{\Gamma }}=\mathrm{log}({M}_{\mathrm{BH}}/{M}_{* })=-1.7$. However, some bias in this observed ratio compared to the intrinsic relation is expected, given the sample selection function (Lauer et al. 2007; Schulze & Wisotzki 2011). In particular, when selecting at the high BH mass end of the relation, the corresponding stellar masses are in the exponentially declining regime of the stellar mass function. This leads to average stellar masses lower than the relation. We estimate the expected bias in Γ following the framework of Schulze & Wisotzki (2011), using the galaxy stellar mass function from Ilbert et al. (2013), and the z = 2 black hole mass function and Eddington ratio distribution function from Schulze et al. (2015). Assuming the SDSS flux limit, and a BH mass selection limit ${M}_{\mathrm{BH}}\gt {10}^{9}$, this predicts a bias of ${\rm{\Delta }}{\rm{\Gamma }}=0.37$ over the Kormendy & Ho (2013) value of ${\rm{\Gamma }}=-2.19$ for black holes with ${M}_{\mathrm{BH}}={10}^{9.5}\,{M}_{\odot }.$ Our quasar hosts are thus consistent with the local relation (${\rm{\Delta }}{\rm{\Gamma }}=0.1$) within uncertainties once selection bias is accounted for, as is the case with previous high-redshift studies (Schulze & Wisotzki 2014).

Zoom In Zoom Out Reset image size

We note in passing that the lack of redshift evolution is consistent with a picture where the scaling relations arise through non-causal means, i.e., galaxies approach the cosmic mean ${\rm{\Gamma }}$ via mergers rather than direct coupling via AGN feedback (Peng 2007; Jahnke & Macciò 2011). However, a lack of evolution does not itself rule out strong AGN feedback, since feedback models can be constructed that predict weak or no evolution.