Host galaxies of high-redshift extremely red and obscured quasars

TL;DR

This study uses Hubble imaging to analyze host galaxies of high-redshift extremely red quasars and type 2 quasar candidates, revealing their luminosities, merger activity, and implications for galaxy evolution models.

Contribution

It provides detailed observations of host galaxy properties and merger activity of high-redshift ERQs and type 2 quasars, challenging existing galaxy formation theories.

Findings

ERQ hosts are very luminous, with some in ongoing mergers.

Type 2 quasar hosts are less luminous and show less merger activity.

Results suggest ERQs may be associated with late-stage mergers, contrasting with some models.

Abstract

We present Hubble Space Telescope 1.4-1.6 micron images of the hosts of ten extremely red quasars (ERQs) and six type 2 quasar candidates at z=2-3. ERQs, whose bolometric luminosities range between 10^47 and 10^48 erg/sec, show spectroscopic signs of powerful ionized winds, whereas type 2 quasar candidates are less luminous and show only mild outflows. After performing careful subtraction of the quasar light, we clearly detect almost all host galaxies. The median rest-frame B-band luminosity of the ERQ hosts in our sample is 10^11.2 L_Sun, or 4 L* at this redshift. Two of the ten hosts of ERQs are in ongoing mergers. The hosts of the type 2 quasar candidates are 0.6 dex less luminous, with 2/6 in likely ongoing mergers. Intriguingly, despite some signs of interaction and presence of low-mass companions, our objects do not show nearly as much major merger activity as do high-redshift radio-loud galaxies and quasars. In the absence of an overt connection to major ongoing gas-rich merger activity, our observations are consistent with a model in which the near-Eddington accretion and strong feedback of ERQs are associated with relatively late stages of mergers resulting in early-type remnants. These results are in some tension with theoretical expectations of galaxy formation models, in which rapid black hole growth occurs within a short time of a major merger. Type 2 quasar candidates are less luminous, so they may instead be powered by internal galactic processes.

Figures

Tables

Table 1. Table 1: Aperture measurements of host galaxies of red and obscured quasars

ID	R.A.	declination	$z$	object	WFC3	host lum.	host lum.	0.3$-$0.6″	0.6$-$1″	1$-$2″
	deg.	deg.		type	filter	$\log \left(\nu L_{\nu }/L_{\odot }\right)$	$\log \left(L_B/L_{\odot }\right)$	flux density	flux density	flux density
SDSSJ0832$+$1615	128.00083	16.25011	2.431	ERQ	F140W	10.90	11.03	8.54 $\pm$ 1.54	3.18 $\pm$ 0.34	0.63 $\pm$ 0.09
SDSSJ0834$+$0159	128.70200	1.98919	2.594	ERQ	F160W	10.89	10.92	8.66 $\pm$ 0.57	2.82 $\pm$ 0.20	0.61 $\pm$ 0.08
SDSSJ0913$+$2344	138.26625	23.74311	2.436	ERQ	F140W	10.92	11.07	17.12 $\pm$ 2.63	2.57 $\pm$ 0.38	0.39 $\pm$ 0.10
SDSSJ1013$+$3427	153.35221	34.45072	2.473	ERQ	F140W	11.02	11.18	12.82 $\pm$ 7.35	4.36 $\pm$ 1.02	1.21 $\pm$ 0.26
SDSSJ1217$+$0234	184.26958	2.57142	2.426	ERQ	F140W	11.41	11.55	21.08 $\pm$ 5.58	7.66 $\pm$ 0.98	1.91 $\pm$ 0.25
SDSSJ1232$+$0912	188.17387	9.20258	2.374	ERQ	F140W	11.20	11.31	13.79 $\pm$ 1.14	5.14 $\pm$ 0.23	2.50 $\pm$ 0.09
SDSSJ1356$+$0730	209.03467	7.50478	2.279	ERQ	F140W	11.17	11.22	16.05 $\pm$ 4.68	2.85 $\pm$ 0.40	0.98 $\pm$ 0.10
SDSSJ1652$+$1728	253.01100	17.48119	2.942	ERQ	F160W	12.03	12.26	53.85 $\pm$ 12.63	15.14 $\pm$ 3.71	14.61 $\pm$ 0.82
SDSSJ2215$-$0056	333.85000	-0.94550	2.493	ERQ	F160W	11.14	11.13	20.45 $\pm$ 1.77	5.39 $\pm$ 0.35	1.03 $\pm$ 0.11
SDSSJ2323$-$0100	350.85904	-1.00919	2.352	ERQ	F140W	11.58	11.68	22.07 $\pm$ 1.16	15.30 $\pm$ 0.25	9.24 $\pm$ 0.08
SDSSJ1001$+$4151	150.39106	41.85220	2.364	T2C	F160W	11.06	11.00	13.11 $\pm$ 2.50	5.42 $\pm$ 0.43	2.01 $\pm$ 0.13
SDSSJ1122$+$3415	170.62645	34.26056	2.398	T2C	F160W	10.58	10.53	6.46 $\pm$ 0.38	1.46 $\pm$ 0.18	0.14 $\pm$ 0.08
SDSSJ1444$-$0013	221.17104	-0.22873	2.548	T2C	F160W	10.36	10.37	3.66 $\pm$ 0.39	1.53 $\pm$ 0.19	-0.06 $\pm$ 0.09
SDSSJ1515$+$1757	228.93337	17.96474	2.402	T2C	F160W	10.90	10.85	15.99 $\pm$ 2.82	3.50 $\pm$ 0.30	0.34 $\pm$ 0.10
SDSSJ2201$+$0012	330.35878	0.20876	2.636	T2C	F160W	11.11	11.17	24.50 $\pm$ 14.53	5.42 $\pm$ 4.92	0.36 $\pm$ 1.21
SDSSJ2233$+$0249	338.45039	2.82579	2.587	T2C	F160W	10.63	10.66	7.80 $\pm$ 0.36	3.12 $\pm$ 0.21	0.39 $\pm$ 0.09

Table 2. Table 2: Results of Sersic fits

ID	$n_s$	$R_{50}$,	$\log$	$\log$	PSF-to
		kpc	$\left(\nu L_{\nu }/L_{\odot }\right)$	$\left(L_B/L_{\odot }\right)$	-host
SDSSJ0832$+$1615	3.2	4.5	11.07	11.21	4.71
SDSSJ0834$+$0159	2.6	5.1	10.66	10.69	9.54
SDSSJ0913$+$2344	1.1	2.6	11.12	11.26	1.29
SDSSJ1013$+$3427			10.70	10.86	$>$10.02
SDSSJ1217$+$0234			11.33	11.46	8.70
SDSSJ1232$+$0912	2.7	6.5	11.14	11.24	3.01
SDSSJ1356$+$0730			10.97	11.02	7.77
SDSSJ1652$+$1728	3.4	3.0	11.76	11.99	2.07
SDSSJ2215$-$0056	2.5	5.2	11.05	11.03	2.14
SDSSJ2323$-$0100			11.64	11.73	0.94
SDSSJ1001$+$4151	3.1	7.3	10.63	10.57	4.80
SDSSJ1122$+$3415	2.8	4.3	10.49	10.44	7.66
SDSSJ1444$-$0013	2.4	4.4	10.23	10.24	7.76
SDSSJ1515$+$1757	2.7	4.9	10.71	10.66	7.24
SDSSJ2201$+$0012			11.32	11.38	$>$ 4.09
SDSSJ2233$+$0249	2.4	5.2	10.60	10.63	4.44

Equations

$$\log I(R)=a_0-(R/b{)}^{1/n_s},$$

$$\log (R_{50}/b)=-0.2963-0.05122n_s+0.2376n_s^2-0.01896n_s^3.$$

$$A=\frac{{\Sigma }_{ij}|I_{ij}-I_{ij,180}|}{2{\Sigma }_{ij}|I_{ij}|}.$$

$$L_{\mathrm{Edd}}=1.3\times 10^{38}\frac{M_{\mathrm{BH}}}{M_{\odot }}erg/sec=1.3\times 10^{47}erg/sec.$$

Full text

††pubyear: 2018††pagerange: Host galaxies of high-redshift extremely red and obscured quasars–Host galaxies of high-redshift extremely red and obscured quasars

1 Introduction

Models for galaxy formation increasingly invoke feedback from active galactic nuclei (AGN) to explain the cessation of star formation in ellipticals and the shape of the galaxy luminosity function (e.g., Silk & Rees 1998; Croton et al. 2006). A popular evolutionary paradigm that naturally incorporates AGN into galaxy evolution models invokes mergers and interactions of gas-rich galaxies to trigger both intense bursts of star formation and black hole growth (e.g., Sanders et al. 1988; Hopkins et al. 2006). In these scenarios, the black holes grow largely in an obscured phase inside the dusty starbursts until a blowout of gas and dust, fueled at least partly by the central AGN, halts the star formation and reveals a visibly luminous quasar in the galactic nucleus. Quasars that are obscured and reddened by dust are valuable for testing this evolution scheme because they should be in the brief blowout phase, close in time to the initial merger/starburst event.

Due to the need for high spatial resolution and high contrast observations in order to detect the faint host galaxies surrounding the bright central quasars, studies of quasar hosts have long been a key science topic for the Hubble Space Telescope (HST). As observational techniques were being developed, pioneering HST studies of quasar hosts naturally concentrated at low redshift ($z<1$, e.g., Bahcall et al. 1997; Kirhakos et al. 1999; Dunlop et al. 2003; Floyd et al. 2004). However, even at low redshifts there remains a significant controversy as to whether quasars are associated with enhanced or decreased star formation in their hosts, whether there are signs of merger activity beyond those expected in similar inactive galaxies, and whether black hole growth is primarily related to stellar mass or star formation rate (Yang et al., 2017). Furthermore, it is unclear how these relationships depend on circumnuclear obscuration, radio loudness or quasar luminosity. For example, Wylezalek et al. (2016) compared type 2 quasar hosts at $z\sim 0.6$ with a stellar-mass matched sample of inactive galaxies, and found that type 2 quasars were hosted by massive galaxies undergoing minor mergers. Canalizo & Stockton (2001) similarly found a connection between quasars in a particular “transition” phase and interactions and star formation in their hosts. Goulding et al. (2018a) found that mergers are particularly important triggers of luminous type 2 quasar activity. But other low-redshift studies, especially of type 1 quasars, have often failed to find a relationship between quasar and enhanced merger activity (Villforth et al., 2017).

Studies of quasar hosts at low redshifts probe the relationship between black holes and their hosts at the epoch when most of the stellar and black hole mass is already in place. In contrast, the epoch $2<z<3$ is particularly interesting for probing models of joint black-hole and host evolution because this is the period when black holes and their hosts were evolving most rapidly (Boyle & Terlevich, 1998; Richards et al., 2006a) and when the quasars must have had the strongest impact on the formation of massive galaxies (Hopkins et al., 2006). But observations of high-redshift active galaxies are even more ambiguous than those at low redshifts, complicated both by cosmological surface brightness dimming and by the relative compactness of high-redshift massive galaxies (van Dokkum et al., 2010).

At $z>2$, the host galaxies of low-luminosity active galaxies in the CANDELS survey (Grogin et al., 2011) appear to be mostly undisturbed and indistinguishable from inactive galaxies of comparable luminosities (Schawinski et al., 2011; Kocevski et al., 2012). These findings echo those of the GOODS survey for active galaxies at $z\sim 1$ (Grogin et al., 2005; Villforth et al., 2012). Perhaps low- and moderate-luminosity sources are not fueled by the merger process and are instead triggered by secular processes within galaxies (Hopkins & Hernquist, 2009), in which case the true tests of merger-driven galaxy and black hole formation scenarios must come from studies of high-luminosity sources.

Until recently, very little was known about the hosts of high-luminosity quasars at $z\sim 2.5$. Attempts to image the hosts of high-redshift luminous quasars were at the limit of detectability (Kukula et al., 2001; Croom et al., 2004; McLeod & Bechtold, 2009), with an occasional aid from gravitational lensing (Peng et al., 2006). Observations of the molecular gas content of powerful quasars using sub-millimeter continuum and molecular line emission has also yielded conflicting results. Some objects appear gas-rich and strongly star-forming (Beelen et al., 2004; Pitchford et al., 2019), while others are devoid of cool gas (Tsai et al., 2015). A slew of recent studies take advantage of the improved resolution and sensitivity of Wide Field Camera 3 (WFC3) on HST, as well as new methods to select obscured quasars, enabling new tests of the merger hypothesis (Glikman et al., 2015; Fan et al., 2016; Hilbert et al., 2016; Mechtley et al., 2016; Farrah et al., 2017). However, disentangling the emission of the bright nuclear source from the surrounding faint and compact host galaxy remains a significant challenge.

In this paper we present the results of an HST study of luminous red and obscured quasars at $z\sim 2.5$, some of which show signatures of extreme quasar feedback. We describe sample selection, ancillary data and HST observations in Section 2. We present host galaxy measurements in Section 3. We discuss our results in comparison with other samples in Section 4, concluding in Section 5. We use a $h=0.7$, $\Omega_{\rm m}=0.3$ and $\Omega_{\Lambda}=0.7$ cosmology throughout this paper. Following long-standing convention, emission lines are identified by their wavelength in air (e.g., [OIII]$\lambda$5007Å), but all wavelength measurements are performed on the vacuum wavelength scale. W1$-$4 refer to magnitudes observed by the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010); we use unWISE flux measurements (Lang et al., 2016). We define the $B$-band luminosity as $\nu L_{\nu}$ in the rest-frame $B$-band and when necessary we apply the conversion for the Sun, $\nu L_{\nu,\odot}$[$B$]$=0.57L_{\odot}$, where $L_{\odot}=3.83\times 10^{33}$ erg/sec is used as a luminosity unit.

2 Sample, observations and data reductions

In this paper, we present results of HST observations of two somewhat different, but related populations of quasars at $2<z<3$. In Section 2.1 we describe selection and ancillary data for a sample of extremely red quasars (ERQs), observed in HST programme GO-14608 (PI: Zakamska). In Section 2.2 we similarly describe selection and ancillary data for a sample of type 2 quasar candidates, observed in HST programme GO-13014 (PI: Strauss). In Section 2.3 we describe our HST observations and data reduction, and in Section 2.4 we describe our modeling of the point spread function (PSF). The targets are listed in Table 1, with spectral energy distributions (SEDs) in Figure 1, separated into the two target categories discussed throughout the paper.

2.1 Sample of extremely red quasars

ERQs were first defined by Ross et al. (2015) who selected spectroscopically confirmed quasars (Pâris et al., 2014, 2017) with $g<22$ mag in the Baryon Oscillation Spectroscopic Survey (Dawson et al., 2013) of the Sloan Digital Sky Survey-III (Eisenstein et al., 2011) based on very red infrared-to-optical colours, as measured from WISE and SDSS photometry. A subset of these objects showed spectroscopic properties not commonly seen in type 1 quasars, such as high ratios of NV$\lambda$1240Å to Ly $\alpha$, high rest equivalent widths (REW) of CIV$\lambda$1550Å and unusual “stubby” (wingless) emission line profile shapes.

Extremely red quasars (ERQs) studied in this paper are drawn from the parent sample by Hamann et al. (2017) who used a combination of selection criteria based on broad-band colours ($i_{\rm AB}-W3_{\rm Vega}>9.8$ mag) and emission-line properties (REW of CIV$>$100Å) to define a sample of 97 ERQs at redshifts $2<z<3.4$. We have been conducting extensive follow-up of these objects across the electromagnetic range. Among the most striking properties of these objects are their extremely high bolometric luminosities, typically $L_{\rm bol}\simeq 10^{47}$ erg/sec, but reaching $10^{48}$ erg/sec in some cases (Hamann et al., 2017), which is an order of magnitude greater than the Eddington limit for a $10^{9}M_{\odot}$ black hole.

Near-infrared (rest-frame optical) spectra of ERQs routinely show forbidden emission lines such as [OIII]$\lambda$5007Å with unprecedented velocity widths, reaching full width at half maximum (FWHM) of 6000 km/sec and showing clear signatures of outflows (Zakamska et al., 2016; Perrotta et al., 2019). Due to the relatively low critical density, this transition is thought to arise on relatively large scales, at least $\sim 100$ pc (Hamann et al., 2011) based on theoretical expectations and out to $\sim 10$ kpc scales (Liu et al., 2013a, b) as seen in low-redshift quasars, although the spatial extents of [OIII] winds in ERQs have not yet been measured.

As a population, ERQs are radio-quiet but not radio-silent (Hwang et al., 2018). Their radio luminosities and [OIII] velocity widths lie on the extreme end of the radio/[OIII] relationship found in low redshift quasars (Zakamska & Greene, 2014). These observations are consistent with the picture in which their radio emission is a by-product of quasar-driven winds, produced when relativistic particles are accelerated in the resulting shocks. Radio data reinforce our hypothesis that ERQs are hosting some of the most extreme galactic winds known (Alexandroff et al., 2016; Hwang et al., 2018).

X-ray observations demonstrate that ERQs are an X-ray luminous, highly obscured population ($N_{\rm H}\simeq 10^{24}$ cm -2; Goulding et al. 2018b) with an X-ray-to-infrared ratio similar to that of other luminous quasars, indicating that most of the infrared emission is likely due to the quasar, not to the star formation in the host galaxy. Although their optical-to-infrared colours are red, their optical (rest-frame UV) continua are relatively bright and flat and their UV emission lines are broad (FWHM of CIV$>$2000 km/sec), inconsistent with those of type 2 quasars (Zakamska et al., 2003). This contradiction between optical and X-ray classifications is resolved by optical spectropolarimetric observations (Alexandroff et al., 2018). These data show that the rest-frame UV continua of ERQs are highly polarized ($\ga 10\%$) and are therefore likely dominated by the light from the obscured quasar scattered into our line of sight by the surrounding medium rather than by direct emission from the nucleus. Furthermore, polarization signatures within the UV emission lines (Ly $\alpha$, CIV) indicate outflow activity on scales of several parsecs. Therefore, outflow activity in ERQs is apparent on a wide range of spatial scales.

For the HST follow-up we prioritized eleven ERQs with then-ongoing follow-up spectroscopic observations covering [OIII]. As near-infrared observations arrived, [OIII] outflow activity in ERQs appeared ubiquitous (Perrotta et al., 2019). Compared to the parent sample of ERQs from Hamann et al. (2017), ERQs targeted for follow-up NIR spectroscopy (and therefore for the HST follow-up) are restricted both in redshift, as [OIII] must fall within the $H$-band or the $K$-band, and in position on the sky for ground-based observability. We have verified that the 11 ERQs selected for the HST study are consistent (in the sense of the Kolmogorov-Smirnov test) with the parent population from Hamann et al. (2017) in all other parameters: infrared luminosity and optical-to-infrared colours, as well as CIV emission line width, equivalent width and kurtosis. The optical-to-infrared SEDs of ERQs selected for the HST observations are shown in Figure 1.

2.2 Sample of type 2 quasar candidates

In the classical geometric unification model of AGN (Antonucci, 1993), obscured (type 2) and unobscured (type 1) objects are distinguished primarily on the basis of the width of the permitted optical emission lines such as H$\alpha$ and H$\beta$ (Khachikian & Weedman, 1974). In this model, the high velocity widths seen in type 1 AGN are appropriate for the gas seen close to the black hole, whereas the low velocity widths seen in type 2 AGN are commensurate with the much shallower potential of the extended host galaxy. A FWHM of 2000 km/sec is typically used to separate type 1 and type 2 quasars (Zakamska et al., 2003; Reyes et al., 2008).

Alexandroff et al. (2013) selected 145 type 2 candidates from the BOSS survey in the redshift range $2<z<4.3$ based on the velocity width of Ly $\alpha$ and CIV. WISE fluxes of these objects indicate bolometric luminosities of up to a few $\times 10^{46}$ erg/sec (Alexandroff et al., 2013). Near-infrared (rest-frame optical) spectroscopy of a subset of these objects (Greene et al., 2014) demonstrates that most of them have a broad H$\alpha$ component and that they therefore have moderate extinction toward the nucleus, $A_{\rm V}\la 2$ mag, although several quasars with type 2 rest-frame optical spectra (i.e., no broad H$\alpha$ or H$\beta$) have also been identified (Alexandroff et al., 2018). Their optical (rest-frame UV) spectra show polarization of a few per cent, suggesting that at least some of the continuum is dominated by scattered light.

[OIII] outflow signatures are present in some of these type 2 quasar candidates (Alexandroff et al., 2018), although they are not nearly as extreme as those seen in ERQs. Given the correlation between line width and radio emission among the radio-quiet quasars (Zakamska & Greene, 2014), it is not surprising that the radio luminosities of high-redshift type 2 quasar candidates are also appreciably lower than those of ERQs (Alexandroff et al., 2016). We conclude that the objects identified by Alexandroff et al. (2013) are in some ways similar to the classical type 2 quasars at low redshifts (Zakamska et al., 2003; Reyes et al., 2008; Yuan et al., 2016), although with a range of extinctions to the nucleus extending to lower values, as appropriate for the selection based on rest-frame UV (instead of rest-frame optical) properties.

As the WFC3/F160W filter fits neatly between the Balmer break and the [OIII]$\lambda$5007Å line for redshifts between 2.4 and 2.6, this was our primary selection criterion of type 2 quasar candidates for HST follow-up. There are 16 objects in this redshift range in the parent sample by Alexandroff et al. (2013). We chose the six objects with the strongest and narrowest emission lines to observe with HST. The FWHM of the CIV line in these objects range from 900 to 1300 km/s, with line luminosities ranging from $8\times 10^{42}$ to $8\times 10^{43}$ erg/sec. Their SEDs from SDSS and WISE photometry are shown in Figure 1. Since they were selected from the SDSS survey without regard for their infrared luminosity, they are not always detected in WISE.

Of the six objects, only two have available near-infrared (rest-frame optical) spectroscopy. One, SDSS J1444$-$0013, shows narrow H$\beta$ and high [OIII]/H$\beta$ ratio characteristic of the narrow-line region, but a strong broad H$\alpha$ component (Greene et al., 2014), so it would be classified as a type 1.9 source based on its rest-frame optical spectrum. The other, SDSS J1515+1757, shows the same type of H$\beta$+[OIII] spectrum (Alexandroff et al., 2018), but has no sign of a broad H$\alpha$ component (Figure 2) and is therefore a spectroscopically confirmed type 2 quasar. In the absence of a spectroscopic type for every object, throughout this paper we continue to refer to this subsample as ‘type 2 quasar candidates’.

2.3 HST images and data reduction

In both HST programmes, we observe each object in the rest-frame ultraviolet with the F814W filter of ACS and rest-frame optical with the F140W or F160W filter of WFC3. At $z\sim 2.5$, the ACS/F814W filter spans rest-frame wavelengths 2000$-$2700Å, an essentially line-free region, thus sampling the near-ultraviolet continuum as a probe both of active star formation and scattered quasar light. The exposure times are one orbit per object for the ERQs and two orbits per object for the type 2 quasar candidates.

The purpose of our near-infrared HST WFC3 observations is to observe starlight from the host galaxy, so for each object we choose the filter to sample longward of the 4000Å break. All type 2 quasar candidates are observed with the F160W filter. In the three lowest redshift objects the [OIII]$\lambda$5007Å emission line falls within the filter and may potentially contaminate the extended emission. For ERQs, the filter choice – F140W or F160W – is adjusted based on the redshift (Table 1), so only in one target (SDSS J2215$-$0056) do we expect the [OIII] line to fall within the filter. In Section 3.4 we discuss the possible contamination of our host observations by emission lines and find it to be negligible. Each target is observed for one orbit with WFC3.

In both programmes, ACS observations are set up with default acs-wfc-dither-box patterns. WFC3 observations of type 2 quasar candidates use both wfc3-ir-dither-box-min and wfc3-ir-dither-blob dithers, while the ERQ observations use only the former. One observation, of ERQ SDSS J0936+1019, suffered from guide star acquisition failure in the WFC3 image (the ACS image is not affected), rendering the image unusable for most of our analyses. Thus the final sample observed with HST consists of 10 ERQs and 6 type 2 quasar candidates, listed in Table 1.

In the following analysis we focus on the WFC3 data to constrain the host galaxy properties, leaving ACS images for future work. We combine the calibrated WFC3 images (flt) of each target into a single drizzled image (drz) using DrizzlePac version 2.1.22. There are 4-6 exposures (flt) per target. In this process cosmic rays and the background are removed and the images are corrected for image distortion and rebinned (“drizzled”) onto a final grid (Fruchter & Hook, 2002). The pixel size of the original flt image for WFC3 IR is 0.13″, but improved resolution can be achieved by taking advantage of the sub-pixel dithers we use in the observations. During final drizzling, we shrink the input pixel (drop size or final_pixfrac) by a factor of 0.8 and subsample the output image with a factor of two finer pixels (0.065″). The final drz images are registered on a coordinate grid with North at the top. We visually inspect the drizzled images and find the drz images to be sharper than the input flt images. In particular, the first Airy ring of the PSF and its substructure are clearly visible in the final images.

2.4 PSF modeling

In Section 3, we present the analysis of the quasar host emission using two different methods – by subtracting the quasar PSF and conducting aperture photometry on the residuals, resulting in 1D surface brightness profiles of the host galaxies, and by jointly fitting PSF+parametric host models to the 2D images. In both cases a PSF model is required. To model the PSF we use Tiny Tim version 7.5, a software package for simulating the HST PSF (Hook et al., 2008; Krist et al., 2011). The fitting is performed on each of the flat-fielded (flt) images that are then combined to make the final drizzled (drz) image. Before the fit, we subtract a constant background level estimated by AstroDrizzle in the original reduction procedure from the flat-fielded image. Then we generate PSF models for the flat-fielded images using Tiny Tim. The PSF model depends on the filter, the focus, the object spectrum, and the position of the object on the detector. We determine the focus of each exposure by fitting the PSF to a reference star in the same field as the quasar, and we determine the spectral type of the reference star from its SDSS colours; the variations in the focus among different images are a negligible source of error. We approximate the spectrum of the quasar by power-law interpolation between the SDSS $z$-band and WISE W1 fluxes.

The PSF model generated has a diameter of 6″ to cover extended spikes and is oversampled by a factor of six to allow sub-pixel shifts. This oversampling factor is chosen as a compromise between being able to capture the structure within the PSF while keeping the computational time manageable. The free parameters are the peak position of the PSF, its normalization, and the width of a symmetric Gaussian kernel that is convolved to the PSF model to simulate pointing jitter. During the fitting process we also apply a charge diffusion kernel provided by Tiny Tim. The best-fit PSF of each exposure is then geometrically corrected and dither-combined with AstroDrizzle in the same way as the images, which accounts for the sub-pixel dithering between constituent exposures. The combined PSF model is then subtracted from the drizzled image (with extension drz).

We find that fits that minimize $\chi^{2}$ often over-subtract the PSF at the nucleus, leaving strong negative residuals in a small number of pixels. This is likely due to the PSF model trying to match the extended features (the host galaxy) by increasing the normalization of the PSF. To reduce the over-subtraction, we penalize the negative residuals by increasing their weights in the $\chi^{2}$ calculation by a factor of 9 relative to the positive residuals.

We use a star in the HST F140W/F160W image of each quasar to evaluate the performance of the PSF subtraction (Figure 3). Within a circular aperture of radius 1″ centered on the star, the PSF-subtracted residual has negligible bias with a standard deviation of 5% of the original flux of the star, whereas beyond 1″, the residuals are $<$10% of the (already small) extended emission of the PSF. Thus our ability to characterize the extended faint emission around bright quasars is limited by our ability to accurately measure the PSF.

3 Host Galaxy Properties

In Section 3.1, we summarize the relationship between the two subsamples and the expectations for the PSF contribution. We detect extended emission in almost all of our 16 targets, as shown in Figures 4 and 5. We conduct two different sets of host galaxy measurements to better understand their uncertainties and the limitations of our observations. In Section 3.2, we discuss our measurements of host fluxes using aperture photometry, in Section 3.3 we discuss model fitting and structural parameters of the host galaxies, and in Section 3.4 we compare the results from the two methods and estimate host luminosities and masses.

3.1 Comparison of type 2 quasar candidates and ERQs

The two subsamples of objects presented in this paper – ERQs and type 2 quasar candidates – were both originally selected from optical and infrared data as candidate luminous obscured quasars. Although the specific selection cuts preclude overlap between the two samples, they show many physical similarities: rest-frame UV continua are likely dominated by scattering (Alexandroff et al., 2018), intrinsic luminosities are high, and rest-frame optical spectra display signs of obscuration commensurate with the selection criteria. Furthermore, both subsamples were observed using a near-identical HST setup allowing direct comparisons. There are also interesting physical differences between these subsamples: ERQs are more luminous and are ubiquitously associated with powerful [OIII] outflows than type 2 quasar candidates, although the origin of the connection between the selection criteria and the [OIII] kinematics is not yet understood (Perrotta et al., 2019).

Although both subsamples are drawn from an obscured population, the HST images have a strong PSF contribution due to the quasars. While we do not perform SED decomposition for our targets (Assef et al., 2015; Farrah et al., 2017), in Figure 1 we illustrate plausible components of the SED and normalize host templates to the median host luminosity in our sample calculated in Section 3.4. Scattered quasar light is likely the source of optical (rest-frame UV) continuum emission, dominating the contribution from the host at these wavelengths. The origin of the near-infrared (rest-frame optical) continuum is less clear. Because broad H$\alpha$ is sometimes detected in these sources (Greene et al., 2014; Zakamska et al., 2016), some direct light from the quasar reaches the observer at these wavelengths, albeit obscured by amounts ranging between $A_{V}<1$ to several mag. Due to the steep rise of the stellar SED at $>0.4$µm, the contrast between the quasar and its host galaxy is most favorable at these wavelengths, which is covered by the WFC3 filters, but PSF subtraction is still required to reveal the quasars’ underlying host galaxies and companions.

3.2 Aperture photometry

After PSF subtraction, for every quasar we calculate aperture photometry in annuli bounded by 0.15″, 0.3″, 0.6″, 1″, 2″, 3″ and 4″ from the center of the PSF. To this end, we sum up the flux within annuli, taking into account the fractional contribution of pixels bisected by apertures. For average surface brightness, we then divide the resulting flux by the area of the annulus. The average surface brightnesses in physical units for the 0.3″-2″annuli are given in Table 1. Surface brightness profiles normalized to the value in the second innermost annulus are shown in Figure 6, left. As the PSF-fitting is performed before the aperture photometry, the results of the fits are independent of the aperture choice.

Next, we use the aperture fluxes of the host galaxies from the PSF-subtracted images to compute the host luminosities. In Figure 6 it is apparent that the surface brightness profiles, over the range where they are well measured, are close enough to $R^{-2}$ that we need to worry about the convergence of the total flux $\int_{0}^{\infty}2\pi R\langle I(R)\rangle{\rm d}R$ both on the low $R$ and on the high $R$ end. Therefore we cannot extrapolate the observed profiles either close to the nucleus or far away from the nucleus where the profiles are not well measured.

To obtain an estimate of the host-galaxy flux, we follow the “minimal flux" approach. We take the PSF-subtracted profile and assume that the surface brightness in the central aperture is the same as in the second aperture (0.15$-$0.3″). For all realistic galaxy profiles the former should be larger than the latter, so by making this assumption we are underestimating the total galaxy flux. Furthermore, the second aperture surface brightness may in turn be under-estimated because of the over-subtraction of the PSF. Then we use the observed data points and integrate the surface brightness profile out to the largest available aperture of 3″, but no further, again potentially underestimating the total flux.

The results are presented in Table 1 in the form of $\nu L_{\nu}$ at the effective rest-wavelength of our WFC3 observations. Furthermore, we bring all luminosity measurements to the same rest-frame Johnson $B$-band111We use a generic Johnson $B$-band filter curve from the Filter Profile Service at the Spanish Virtual Observatory http://svo.cab.inta-csic.es/ by applying k-corrections using the Sb templates from Polletta et al. (2007) shown in Figure 1 and hereafter report $L_{B}\equiv\nu L_{\nu}$ at rest-frame 4400Å. In the wavelength range 0.4$-$0.7µm, Sb and Ell templates produce similar results, whereas the M82 template, which has more dust obscuration, is appreciably redder in the optical. Because our observations are close to rest-frame $B$-band, k-corrections are in general smaller than other uncertainties – out of the objects in our sample, two thirds have k-corrections under 0.1 dex and the remaining third under 0.2 dex. The median measured value is $L_{B}=10^{11.1}L_{\odot}$, and the mean and the sample standard deviation of $\log(L_{B}/L_{\odot})$ are $11.1\pm 0.4$.

Finally, we can estimate the effective sizes and measures of concentration by comparing the observed profiles with Sersic profiles. A circularly symmetric Sersic profile is a three-parameter function, with surface brightness given by

[TABLE]

where $a_{0}$ is the overall normalization and $n_{s}$ is the Sersic index: $n_{s}=1$ for an exponential profile (Patterson, 1940) and $n_{s}=4$ for a de Vaucouleurs (1953) profile.

The parameter $b$ is proportional to the physical scale of the galaxy, which historically is defined differently for different Sersic indices. One physically meaningful non-parametric value is effective radius $R_{50}$ – the radius containing half of the total galaxy light. The relation between $R_{50}$ and $b$ can be calculated numerically for every $n_{s}$, and in the range $n_{s}=0.5-4$ it is well approximated by a third-order polynomial:

[TABLE]

For $n_{s}=4$ (de Vaucouleurs), this equation yields $(R_{50}/b)^{1/4}=3.33$, which corresponds to the classical de Vaucouleurs’ parametrization.

Thus, a circularly symmetric Sersic profile can be completely described by the normalization, index $n_{s}$ and half-light size $R_{50}$, whereas flux-normalized Sersic profiles form a two-parameter family. In the right panel of Figure 6 we show model surface brightness distributions obtained from Sersic profiles using the same method as in the observations: i.e., convolving with the PSF, integrating within annuli and then dividing by their areas. We probe a plausible range of Sersic indices ($n_{s}=1,2,4$) and effective radii $R_{50}$ of 0.25, 0.5, 1.0 and 2.0″. Qualitatively, most of the observed surface brightness profiles appear to be on scales similar to those of Sersic profiles with $R_{50}=0.25-1\arcsec\simeq 2.1-8.4$ kpc, which is unsurprising considering that that most massive galaxies at $z>2$ have $R_{50}<3$ kpc (van Dokkum et al., 2010).

3.3 Sersic fits and structural parameters

To obtain another measure of the host flux and determine the structural parameters of the host galaxies in our sample, we perform a joint PSF + two-dimensional Sersic fit, which should alleviate the problem of over-subtraction in our previous PSF-only fits. There are many challenges to this analysis. First, as the hosts may be compact, some of their light can be hidden within the central PSF, so we expect some degeneracy between the contributions of the host galaxy and of the quasar. This is a generic difficulty of modeling quasar hosts in imaging data, including at low redshifts (Zakamska et al., 2006). Second, as discussed in Section 2.4, the accuracy with which the PSF can be modelled is $\sim 5\%$ in the center and $\sim 10\%$ in the outer parts, limiting the quality of the fits for our sources, as the PSF is larger than the flux of the host galaxy by up to an order of magnitude (Figure 1).

While fitting the flt images, we keep all the parameters of the PSF and Sersic models free. For the PSF model the parameters include its position, amplitude, a Gaussian smoothing width to simulate jitter, and a background level. The two-dimensional Sersic model which represents the host is described by 7 parameters: coordinates of the central position, Sersic index, normalization, half-light radius $R_{50}$, ellipticity and position angle. The equations connecting these parameters to total fluxes and other standard parametrizations of the Sersic profiles are listed in Section 3.2 and in Zakamska et al. (2006). To avoid over-fitting, the Sersic indices are constrained to the range $0-4$ and the ellipticities to the range $0-0.4$. During the fitting, the two-dimensional Sersic model is convolved with a two-dimensional Gaussian best fit to the PSF model.

The best-fit PSF models of the flt images are drizzle-combined to form the PSF model of the drz image. This PSF model is then fitted to the drz image with the amplitude as the only free parameter, together with a free Sersic model. The central residuals can be quite large, higher than the galaxy signals, due to both the Poisson noise and the mismatch between the PSF template and the data, so they contribute overwhelmingly to the $\chi^{2}$. To avoid placing all statistical weight of the fits on the region where the host signal contributes the least, the central pixels that are brighter than 20% of the peak intensity are masked at this stage. Companions, foreground stars, and galaxies are masked so that they do not affect the orientation and the ellipticity of the Sersic profile. Because of the higher signal-to-noise ratio and the better resolution of the drz-based images, we use drz-based fits with these masks in what follows. Some of these fits are shown in Figure 7, and the fitting parameters are presented in Table 2.

In most cases the joint PSF+Sersic model fits the data reasonably well, leaving only small residuals that may be due to the asymmetry of the host galaxy. There are a few cases where the diffraction spikes of the PSF are not entirely subtracted by the best-fit model, including J0834+0159, J1013+3427, J1217+0234, J1652+1728, and J2201+0012. If we increase the amplitude of the PSF until the spikes are subtracted, the core of the PSF model becomes brighter than the original image. These few targets also have the highest PSF intensity among the sample, suggesting that the core of these images have a nonlinear or saturated response. Residual diffraction spikes and rings would increase the fluxes on the outskirts of the galaxy, thus lowering the measured Sersic index, increasing the radius and leading to additional uncertainties in the flux.

There are multiple sources of systematic uncertainties in our fitting results, including imperfections of the PSF model (e.g., possible saturation issues) and deviations between the real surface brightness profiles and the model assumptions – for example, that the host galaxy is well represented by a Sersic profile. Such systematic uncertainties dominate the statistical error. We test the sensitivity of the outputs to the assumptions about the PSF in a variety of ways by modifying the fitting procedure or the masking prescriptions. For example, we find that flt-based fluxes are $\sim 0.2$ dex brighter than the drz-based ones, which is likely due to the lower effective resolution of the flt-based images and the resulting higher degeneracy between the nucleus and the host. The standard deviation of the difference between flt-based and drz-based fluxes is 0.25 dex, whereas the variations among the fitting procedures we have tried result in standard deviations $\la 0.2$ dex. We adopt the drz-based fluxes, and use 0.25 dex as a conservative estimate of the Sersic flux error. In Table 2 we list both the measured monochromatic luminosities and those k-corrected to the rest-frame $B$-band, finding for the entire sample $\log(L_{B}/L_{\odot})=11.0,11.0\pm 0.5$.

Furthermore, we find that the aperture fluxes and the Sersic fluxes show no significant offset relative to one another, with aperture-based fluxes being are on average higher by only 0.1 dex than the Sersic-based fluxes. The standard deviation of the differences between the two sets of fluxes is 0.18 dex, which is within our estimated error on the Sersic fluxes. The two most discrepant fluxes are found in hosts with companions, which are naturally not well modeled by a single Sersic profile. Excluding the five objects with companions named in Sec. 3.5 reduces the standard deviation of the difference to 0.14 dex. This comparison of fluxes calculated by two completely different methods suggests that despite the challenging PSF subtraction, the extended host flux measurements are quite robust. However, if the host galaxy has a compact circumnuclear starburst, it would be absorbed into our PSF model and not reported as part of the host flux with either the aperture or Sersic method, and we cannot make a realistic estimate of uncertainty associated with this possibility. Emerging integral field unit spectroscopic studies of quasar hosts may be able to circumvent this problem by taking advantage of spectral differences between the quasar and its host (Rupke et al., 2017).

We display the distributions of our best-fit Sersic indices $n_{s}$ and half-light radii $R_{50}$ in Figure 8 and list them in Table 2. We expect significantly higher relative uncertainties in the structural parameters than in the flux because $n_{s}$ and $R_{50}$ are somewhat degenerate, with both high $n_{s}$ and small $R_{50}$ resulting in a compact light distribution. To quantify this uncertainty, we simulate two-dimensional Sersic profiles with a fixed total flux outside of the inner 0.15″ (mimicking our PSF contamination). We then find a family of Sersic models with surface brightness at 0.5″, 1″ and 2″ within 10% of the fiducial model with $n_{s}=2.5$ and $R_{50}=4.5$ kpc. Within this family, whose surface brightness profiles would look essentially indistinguishable to an observer, $n_{s}$ ranges between 1.8 and 4 – i.e., it is very poorly constrained. The range of $R_{50}$ is narrower, between 3.5 kpc and 5 kpc (relative error of $\sim 0.1$ dex), due to its stronger sensitivity to the flux constraints outside 0.15″.

An additional source of error in structural parameters is due to the fitting itself, which we estimate by comparing results from different fitting procedures (e.g., with and without masks), from different flt images, and from drz vs flt images. The resulting fitting uncertainty for each source (estimated as the sample standard deviation from a range of fits) is typically 0.5$-$0.6 for $n_{S}$ and 0.18 dex for $R_{50}$. As with fluxes, if there is a bright centrally concentrated component (high $n_{s}$ or small $R_{50}$), the contribution from the central part of the galaxy is degenerate with the quasar PSF, so all information about such component is lost due to the presence of the quasar, and the systematic uncertainties due to this effect cannot be estimated. Given all these considerations, we adopt a nominal error of 1.0 on each $n_{s}$ measurement and 0.2 dex on each $R_{50}$ measurement except for the 5 most compact objects ($R_{50}\la 4.5$ kpc) where the uncertainty is 0.3 dex.

A Kolmogorov-Smirnov test shows that the distributions of these structural measurements for ERQ and type 2 quasar candidate hosts are consistent with being drawn from the same population. The fit values of $R_{50}=4.9,4.8\pm 1.3$ kpc or $0.6,0.6\pm 0.15\arcsec$ (median, mean and sample standard deviation) are qualitatively consistent with our expectations from Figure 6. When compared to local galaxies, the Sersic indices $n_{s}=2.7,2.6\pm 0.6$ are intermediate between disk-like and bulge-like, but when compared to distant massive galaxies, these values are typical (van Dokkum et al., 2010). However, while our estimated structural parameters appear to be similar to those of high-redshift massive galaxies, the uncertainties on these values are too large at this point to warrant a more detailed comparison with the high-redshift populations studied with much deeper data and without the difficulties of the central source contamination. Furthermore, because the individual errors on $n_{s}$ and $R_{50}$ are comparable to the standard deviation within the population, the measurements are not yet sensitive to any differences between the structural parameters of ERQ and type 2 quasar candidate hosts.

The three objects with measured effective radii in excess of 10 kpc are SDSS J1217+0234, SDSS J1356+0730 and SDSS J2323$-$0100, all of them ERQs. SDSS J2323$-$0100 is an ongoing merger, where a single-Sersic fit is inadequate. In the other two objects the host galaxy is clearly detected in the PSF-subtracted images, but is not so extended as to warrant such large effective radii. It is likely that either a Sersic profile is an inadequate representation of the galaxy, or the PSF model is not very accurate, and we consider these joint fits to be unsuccessful and omit their structural parameters from Table 2.

3.4 Host luminosity and stellar mass

In what follows we use $L_{B}$ luminosities from Sersic fits, but none of the results are sensitive to which set of luminosities we choose to use. The median luminosity of all host galaxies in our sample is $L_{B}=10^{11.0}L_{\odot}$. The mass-to-light ratio of stellar populations in the $B$-band is a function of metallicity, age and star formation history. For stellar populations with age $0.3-2$ Gyr (the upper limit being the maximum age at the redshift of our sample) and solar metallicity $Z_{\odot}$ characteristic of our targets (Hamann et al., 2017), the mass-to-light ratio ranges between 0.4 and 4 $M_{\odot}/L_{\odot}$ (Maraston, 2005), where we have corrected from the $B$-band luminosity used by Maraston (2005) to the bolometric luminosity $L_{\odot}$ used as a luminosity unit in this paper. Therefore, the median corresponding stellar mass is $10^{10.4-11.4}M_{\odot}$, depending on the age of the galaxies in our sample. A higher mass-to-light ratio is possible if the galaxies are very dusty, and an 0.3 dex lower mass-to-light ratio would be expected for $Z_{\odot}/20$ – e.g., if the host galaxy has a much lower metallicity than the nucleus (Maraston, 2005).

The measured galaxy fluxes imply that the extended narrow-line region contamination to our imaging (Section 2.3) is minor. The total [OIII] luminosity rarely exceeds $10^{10}L_{\odot}$ (Greene et al., 2014; Zakamska et al., 2016), which is less than 10% of the host luminosity falling within the broad filter. Furthermore, only a small fraction of the line luminosity is expected on scales of a few kpc, as its spatially unresolved, centrally-concentrated component is subtracted using our PSF subtraction procedure. Thus the extended line luminosity makes a negligible contribution to our measured extended host fluxes. In order to detect extended line emission, one would need to use integral field spectroscopy (Liu et al., 2013a) or narrow-band imaging.

Despite limited statistics we find an appreciable luminosity difference between the hosts of ERQs and those of type 2 quasar candidates in our sample (Figure 8, right). The median, the mean and the standard deviation of the log host luminosity for the ERQ sample are $\log(L_{B}/L_{\odot})=11.2,11.3\pm 0.4$, whereas for the type 2 quasar candidate sample these values are $\log(L_{B}/L_{\odot})=10.6,10.7\pm 0.4$. The Kolmogorov-Smirnov test shows that the two distributions are statistically inconsistent; the probability that the two host luminosity distributions are drawn from the same underlying one is 0.4%. Since each distribution is likely a convolution of the intrinsic distribution with the measurement errors of $\sim 0.25$ dex, accounting for errors would make the apparent difference even more statistically significant. Thus we find evidence that ERQ hosts are either more massive or more strongly star-forming than the hosts of type 2 quasar candidates, resulting in a greater $B$-band luminosity.

For comparison, the $B$-band luminosity $L^{*}$ of the break of the galaxy luminosity function at $z=2.5$ is at $m_{AB}=-21.74$ mag or $10^{10.6}L_{\odot}$ (Giallongo et al., 2005). Therefore, type 2 quasar candidates’ hosts have luminosities right around $L^{*}$ at their redshift, but ERQ hosts are on average four times more luminous. In Section 4.1 we demonstrate that there is no evidence for scattered light contamination of the extended emission in the rest-frame $B$-band. Another possible explanation for the difference in host luminosities is that type 2 hosts might suffer from a greater level of galactic-scale extinction than do ERQs. While we cannot conclusively rule out this possibility, we see no evidence for any differences in structural parameters or morphology between the two subsamples. Thus in what follows we take the difference between ERQ and type 2 host luminosities at face value.

The host-to-PSF flux ratio is $0.4,0.3\pm 0.3$ for the ten ERQs (median, mean, sample standard deviation) and $0.18,0.17\pm 0.05$ for the six type 2 candidates (the large standard deviation of the ratio for the ERQs is partly because SDSS J2323-0100 has a much larger host-to-PSF ratio than the rest of the sample). The consequence is that the nuclei of the type 2 quasar candidates in our sample are nearly an order of magnitude less luminous in the rest-frame optical than the nuclei of ERQs. This is not surprising: the majority of ERQs show evidence for a broad-line region in their optical spectra (Zakamska et al., 2016; Perrotta et al., 2019) and therefore have a type 1 optical classification. In contrast, type 2 quasar candidates typically have a type 2 or a type 1.9 classification (Section 2.2) and therefore a smaller contribution of the nuclear light to the rest-frame optical spectra.

3.5 Quantifying signs of merger activity

Models of co-evolution of black holes and their hosts often appeal to mergers as a leading agent in triggering the active nucleus and the formation of the spheroid (Hopkins et al., 2006). Some models predict that nearly 100% of quasars with $L_{\rm bol}>10^{47}$ erg/sec reside in galaxies that are merging, interacting or irregular (Treister et al., 2012; Hickox et al., 2014). However, there is no universal agreement on how to quantify such activity in quasar hosts, and identifying major merger activity in quasar hosts at high redshifts remains an especially difficult problem.

At low redshifts, one of the most identifiable hallmarks of galactic interactions is the presence of tidal tails, but these are hard to see at $z=2.5$ due to cosmological surface brightness dimming. Furthermore, the density of sources in WFC3 images is quite high, and two galaxies that are close together on the plane of the sky can be a chance superposition. Finally, the disturbed morphology of the host may be a sign of clumpy star formation characteristic of high-redshift galaxies rather than a recent merger (Agertz et al., 2009), especially in the rest-frame $B$-band probed by our WFC3 observations. Now that several studies of hosts of extremely luminous quasars are available in the literature, our primary interest is to determine whether the wide variations in the reported merger fractions are related to any quasar properties. The comparison among different studies is complicated by the use of widely different methods – often just visual classification without any PSF subtraction (e.g., Fan et al. 2016) – in the literature.

Keeping this in mind, we classify a target quasar as an ongoing major merger if it has a galaxy within 2.5″ with a flux that is greater than 1/4 of the estimated flux of the quasar host galaxy (Goulding et al., 2018a; McAlpine et al., 2018). While this definition likely provides a lower limit on major merger activity, it has the advantage of being nearly redshift-independent, enabling us to compare studies in the range $1<z<3$. Indeed, it does not rely on diffuse features, such as tidal tails or shells. For comparison, over 90% of local ultraluminous infrared galaxies (ULIRGs) with a double-nucleus morphology have a flux ratio above 1/4 (Kim et al., 2002). Furthermore, the angular diameter distance is within 5% of the $z=2.5$ value in this range. Finally, by requiring a high flux ratio between the host and the companion and a small separation distance, we reduce the probability of counting chance superpositions and small companions without dynamical influence on the host (McAlpine et al., 2018). The major drawback of this definition is that a late-stage merger without two separate galactic nuclei would not be identified as such. For example, only 42% of the local ULIRGs have a double-nucleus morphology we would recognize in our observations, with the rest being single late-stage remnants (Kim et al., 2002).

While there are many faint objects in our images, only three targets out of 16 have a companion within 2.5″ with a flux above 1/4 of our Sersic-fitted flux: two are type 2 quasar candidates (SDSS J1001+4151 and SDSS J1444$-$0013) and one is an ERQ (SDSS J2323$-$0100). ERQ SDSS J1217+0234 has a companion (itself a clumpy, disturbed galaxy) with a flux of 0.24$\times$ the flux of the host galaxy, just under the 1:4 cutoff, and could be a major merger within the uncertainty of our flux measurements. For each image, we calculate the probability of a chance superposition within 2.5″ to be 0.11 using WFC3 number counts by Guo et al. (2013) and the range of fluxes between those of our targets and our cutoff for the companion (1/4). We further use WFC3/ACS flux ratios and redshifted Polletta et al. (2008) templates to confirm that the identified companions are at $z>1$ (the ACS images lack the depth for better redshift constraints), which decreases the probability for a chance superposition by about a factor of 2.

In addition to the major merger candidates selected due to the presence of a companion, we find an additional candidate – ERQ SDSS J1652+1728. The analysis of this target is complicated by the presence of a bright compact source just North of our object (Figure 9); from the combination of optical colours from the SDSS, our derived HST near-infrared colours and Keck spectroscopy (Alexandroff et al., 2018) which shows a featureless spectrum, we conclude that this object is most likely a foreground star. Excluding this chance superposition from consideration, we find that SDSS J1652+1728 has multiple faint companions within 2″ – which could well be giant clumps of star formation characteristic of high-redshift star-forming galaxies (Elmegreen et al., 2007; Förster Schreiber et al., 2009) – and a pronounced tidal tail extending over tens of kpc, marked with an arrow in Figure 9. Therefore, we flag this object as a major merger candidate even though it does not meet our formal companion-based definition. Thus our final count of ongoing major mergers among ERQs is $2\pm 1$ out of 10 objects, and it is 2/6 for the type 2 quasar candidates. There are no other strongly disturbed sources that we would be tempted to classify as major mergers, although a late-stage single-nucleus merger remnant can be easily missed by our classification.

In a late-stage merger, a separate companion would not necessarily be identifiable and tidal features fade, but the host galaxy could still show morphological disturbances. Various non-parametric quantitative measures of morphology exist, such as $M_{20}$, Gini coefficient and asymmetry (Lotz et al., 2004) defined as

[TABLE]

Here $I_{ij}$ is the original 2D image, while $I_{ij,180}$ is the image rotated by 180o about the central point – in our case, the quasar. The key complication in using these parameters for our sample is the PSF subtraction. In particular, both $M_{20}$ and Gini coefficients rely upon the relative distributions of bright and faint galaxy light. Since the bright center of the host galaxy is blocked by the PSF, these measures are not suitable for our sample where the PSF subtraction uncertainties dominate these regions. The asymmetry is more promising, because if the PSF is symmetric, its residuals would not contribute to the numerator of equation (3).

However, we find that the asymmetry measurements vary by over a factor of 2 depending on the exact method of PSF subtraction and to the range or radii over which they are calculated. The same difficulty was also pointed out by Glikman et al. (2015). We conclude that non-parametric measurements of morphology are not yet well tailored to galaxies with a contaminating central source. We use asymmetries with extreme caution and only for comparison with Glikman et al. (2015) when we follow the same measurement procedure as they do. This involves computing asymmetry within 1.5$\times$ the Petrosian (1976) radius calculated from our Sersic fits, while masking a 1″$\times$1″region around the nucleus and not masking the companions, resulting in median, mean and sample standard deviation of asymmetries of 0.58, 0.53$\pm$0.15. While ERQs and type 2 quasar candidates are statistically indistinguishable in this parameter as per the Kolmogorov-Smirnov test, this may not be a meaningful comparison given the restrictions imposed on our measurement procedure and large individual errors. We discuss the implications of our morphological measurements in the next section.

4 Discussion

With our host galaxy measurements in hand, we discuss the implications of our results and compare the hosts of our targets with several other samples. In Section 4.1 we determine whether light scattered off the interstellar medium of the host galaxy can affect our measurements. In Section 4.2 we discuss the black hole vs host galaxy relationship of our sample and in Section 4.3 we discuss the constraints on the star formation rates. In Section 4.4 we summarize several comparison samples of powerful active nuclei at similar redshifts with host galaxy measurements. We compare the hosts of these objects with ours in Sections 4.5 and 4.6.

4.1 Effects of scattered light

Even when the direct view to the quasar is obscured – either by the host galaxy or by circumnuclear dust – the emission from the quasar can escape along other directions, scatter off the surrounding interstellar matter and reach the observer (Antonucci & Miller, 1985; Antonucci, 1993). Alexandroff et al. (2018) present Keck spectropolarimetric observations of ERQs and type 2 candidates, four of which overlap with our HST observations. Three of these four – SDSS J1232+0912, SDSS J1515+1757 and SDSS J1652+1728 – show high continuum polarization in the optical (rest-frame UV), between 10 and 15%, while the remaining source SDSS J2215$-$0056 is too weak in the continuum to reliably measure polarization. While these values are not at the theoretical maximum for polarization of scattered light, they are at the high end of polarization for the obscured quasar population (Hines & Wills, 1993; Smith et al., 2002; Zakamska et al., 2005; Greene et al., 2012). Therefore, if these three sources are representative of the entire HST sample, then in principle all of the optical (rest-frame UV) emission can be due to scattered light.

As the wavelength increases, the scattered light contribution decreases while the host galaxy becomes more important, so it is not clear what fraction of the flux and on what scales the scattered light contributes to our WFC3 images. If the extended emission outside of the PSF was dominated by scattered light, we would expect a 90 degree offset between the position angles of the optical emission and the position angle of the E-vector of polarization. Instead, for the first three objects we find differences of $\sim 20$, $\sim 50$, and $\sim 30$ degrees. In the last target neither position angle is well determined, but our best attempt yields $\sim 40$ degrees. The clearest measurement is in SDSS J1515+1757, with Figure 7 demonstrating that the position angle of the major axis of the fitted galaxy is 30 degrees E of N, whereas Figure A2 of Alexandroff et al. (2018) shows that the reddest polarization position angle measurement is 80 degrees E of N, leading to the offset between the two of 50$\pm$10 degrees.

Thus none of the objects for which we have both HST and polarimetric data show the 90 degree offset expected if the extended emission was dominated by scattered light. This is somewhat surprising in light of the results for lower-redshift obscured quasars (Zakamska et al., 2006; Obied et al., 2016), in which the polarization position angle and the extended blue emission are often accurately orthogonal to each other. Therefore in contrast to the low-redshift obscured quasars, the scattering must predominantly occur in a high-density medium on scales smaller than can be probed with our WFC3 observations ($\la 1$ kpc). In this case the scattered light would be part of the PSF component of our analysis. This conclusion is consistent with the estimates of the spatial scales of scattering presented by Alexandroff et al. (2018), who found that the rest-frame UV emission lines must originate on scales that are similar to the scales of scattering, $\la$ a few tens of pc. Thus, even though the scattered light constitutes a major component of the observed SEDs of our sources, it does not appear that it has much impact on the host galaxy properties derived from the WFC3 data. Extended scattered light can be further probed by a careful analysis of the ACS images where the contrast between the scattered light and the host galaxy is more favorable.

4.2 Black hole masses and the nature of our targets

What kinds of black holes power the extremely luminous quasars at these redshifts? While there exist established scaling relationships between the host stellar mass and black hole mass in nearby galaxies (Kormendy & Richstone, 1995), comparing them with our objects presents major difficulties. First, there is a potential evolution in the black hole vs galaxy relationship between the epoch corresponding to our targets and the present day. Second, we are observing in the rest-frame $B$-band, which is strongly affected by on-going and recent star formation and may be a poor probe of stellar mass. Third, given the uncertainties in our measured Sersic indices and the fact that their values are intermediate between disc-like and bulge-like, it is not clear whether we are detecting the progenitors of massive bulges.

With these caveats, we concentrate on our ERQ subsample, where the mid-infrared luminosities are well measured, giving us a good handle on the bolometric luminosity. To avoid systematic errors from conversions between different bands and stellar masses, we use local galaxies with $B-$band luminosities from Savorgnan et al. (2013), but they have no objects in their sample with $L_{B}>10^{11.2}L_{\odot}$. Therefore, the best available comparison subsample of local galaxies from Savorgnan et al. (2013) has $L_{B}=10^{10.8-11.2}L_{\odot}$, with a median luminosity 0.3 dex less luminous than ours. The median, mean and standard deviation of the dynamically measured log black hole masses within the comparison subsample are $\log(M/M_{\odot})=8.9,8.7\pm 0.9$.

With these crude estimates of black hole mass in local galaxies with $L_{B}$ luminosities comparable to those of ERQs, we confront three different possibilities. The first is that the black hole masses of ERQs are also $\sim 10^{9}M_{\odot}$, in agreement with the local black hole vs host relationships. The rest-frame 5µm luminosities of ERQs of $10^{46}-10^{47}$ erg/sec (Figure 1) correspond to bolometric luminosities of $10^{47}-10^{48}$ erg/sec (Richards et al., 2006b). At a $10^{9}M_{\odot}$ black hole mass, ERQs would have to be up to 10 times brighter than the Eddington luminosity:

[TABLE]

One possible argument against such high inferred Eddington ratios is that we might be over-estimating bolometric luminosities by using the bolometric corrections characteristic of type 1 quasars from Richards et al. (2006b). However, this is unlikely because in the X-ray vs infrared luminosity space ERQs are close to type 1 quasars (Goulding et al., 2018b). In fact, for an obscured population it is more plausible that the bolometric corrections are under-estimated rather than over-estimated, and somewhat higher bolometric corrections would still be consistent with the X-ray data (Goulding et al., 2018b).

The structure of the outflows in the near-Eddington and super-Eddington regime is sensitive to the complex radiative transfer in the optically thick regime. Strong radiation pressure may lead to powerful outflows (Shakura & Sunyaev, 1973), or photons can be trapped in the matter and advected into the black hole (Begelman & Meier, 1982). This regime has only recently been probed by numerical simulations (Sądowski & Narayan, 2016; Dai et al., 2018), which demonstrate that at high accretion rates, wide-angle outflows ubiquitously form and carry an appreciable fraction of the accretion power. The extreme wind activity seen in our targets on multiple scales (Zakamska et al., 2016; Alexandroff et al., 2018; Perrotta et al., 2019) is qualitatively consistent with super-Eddington accretion, so our multi-wavelength observations of ERQs support this scenario, and other populations of candidate super-Eddington accretors also show powerful outflow activity (Done et al., 2007; Ni et al., 2018). Furthermore, simulations confirm earlier suggestions that due to photon advection, near- and super-Eddington accretion tends to be radiatively inefficient (Abramowicz et al., 1988). Therefore, to produce observed ERQ luminosities of a few times the Eddington limit, an even higher ratio of the accretion rate to $\dot{M}_{\rm Edd}$ may be required.

Another possibility is that the black holes in ERQs are appreciably more massive, with $M_{\rm BH}\sim 10^{10}M_{\odot}$, bringing their luminosities just below the Eddington limit. This might occur if the local black hole vs host relationship still holds for ERQs, but the $B$-band luminosities of our sample are under-estimated by a factor of ten (e.g., due to strong host galaxy obscuration). Alternatively, the $B$-band host luminosities are correct, but the black holes are ten times more massive than suggested by the local black hole vs host relationship. Both these options require extremely high black hole masses, at the upper limit of those observed in the local universe (McConnell et al., 2011; Thomas et al., 2016). This would imply that in our target ERQs, the observed quasar activity is the last episode of appreciable black hole growth. With a better understanding of the duty cycle of ERQs, the different hypotheses about the masses of their supermassive black holes can be tested by measuring the volume density of ERQs and comparing it to the local volume density of black holes of different masses.

4.3 HST constraints on star formation

Here we again concentrate on ERQs whose SEDs are better measured than those of type 2 quasar candidates. The star formation rates of ERQ host galaxies are still not well-constrained with the existing data. If their radio luminosity were due to star formation in the host galaxies, then the observed median radio luminosity implies an extremely high star formation rate of 2700$M_{\odot}$/year (Ross et al., 2015; Alexandroff et al., 2016; Hamann et al., 2017). Another potential argument in favor of extremely high star formation rates is supplied by the model of co-evolution of active galaxies and their hosts by Hickox et al. (2014). These authors find that for $L_{\rm bol}=10^{47-48}$ erg/sec quasars at $z=2$ the expected infrared luminosity of star formation is in the range $10^{46.2-47}$ erg/sec, nominally corresponding to star formation rates of $700-4500$ $M_{\odot}$/year (Bell, 2003).

Such extremely high star formation rates are well above those seen in sub-millimeter galaxies (Blain et al., 2002; Chapman et al., 2005), but have been seen in the most extreme high-redshift starbursts (Rowan-Robinson et al., 2018) and may not necessarily be out of place in the hosts of ERQs, which are after all some of the most rapidly accreting active nuclei in the universe. Nonetheless, we previously argued that such high star formation rates are unlikely because the radio luminosities of ERQs are dominated by the active nuclei, not by the host galaxies (Hwang et al., 2018). Furthermore, the co-evolution models are not yet well probed at $L_{\rm bol}>10^{47}$ erg/sec, and the existing observations do not support the predicted extremely high star formation rates at the high luminosity end (Hickox et al., 2014).

Thus the star formation rates in ERQ hosts remain controversial. We now consider whether our HST data provide any additional constraints on the star formation rates. For the nominal star formation rate of 2700$M_{\odot}$/year necessary to account for the radio emission, in the rest-frame $U$-band we expect a luminosity of $L_{U}\equiv\nu L_{\nu}[U]=6\times 10^{12}L_{\odot}$ due to hot young stars (Moustakas et al., 2006). For an O5 stellar template (Pickles, 1998), we find $L_{B}/L_{U}=0.6$, and this ratio is higher for all later stellar types. Therefore, the median $B$-band luminosity at the star formation rates necessary to explain the observed radio emission and to bring our objects in agreement with models by Hickox et al. (2014) is 1.4 dex higher than observed in our HST data. This is yet another argument against extremely high star formation rates in ERQs. However, much of star formation could proceed in an obscured fashion, as seen in other obscured quasars at the same redshifts (Wethers et al., 2018); therefore, we still need alternative measures of star formation in our sources, e.g., from the far-infrared or submm continuum, to settle this issue.

4.4 Introducing comparison samples

Despite the strong contribution of the quasar PSF to the observed emission, we successfully detect almost all host galaxies of our target quasars. The ability to routinely detect host galaxies of powerful quasars at $z>2$ is a major success of WFC3 sensitivity, image quality and improved analysis tools, and there are not yet many samples available for comparison with ours. Here we focus on the samples that are most directly comparable in terms of quasar luminosity, redshift and observing mode, but they span a range of obscuration and other properties.

Mechtley et al. (2016) analyze host galaxies of 19 radio-quiet type 1 (unobscured) quasars at $z=1.9-2.1$ selected based on their luminosity and black-hole mass. The median optical luminosity of the quasars in this sample is $L_{\rm 3000\AA}=10^{46.4}$ erg/sec or $L_{\rm bol}\simeq 10^{47.1}$ erg/sec, using bolometric corrections by Richards et al. (2006b). These objects are observed for one orbit with the F160W filter of WFC3, and among other analyses Mechtley et al. (2016) perform PSF+Sersic decomposition and merger classification of their samples.

Fan et al. (2016) and Farrah et al. (2017) analyze one-orbit F160W or F110W images of hot dust-obscured galaxies (HotDOGs) at $z=1.8-3.9$ which are selected based on their red WISE colours (Eisenhardt et al., 2012; Tsai et al., 2015; Assef et al., 2015). X-ray observations indicate that these objects host powerful obscured quasars (Assef et al., 2016; Ricci et al., 2017; Vito et al., 2018) which likely dominate the bolometric output (Farrah et al., 2017). HotDOGs have similar luminosities and shapes of spectral energy distributions (SEDs) to ERQs (Goulding et al., 2018b), but HotDOGs likely have higher obscuration levels as their SED rise is at longer wavelengths than in ERQs (Hamann et al., 2017). Furthermore, in the admittedly limited sample of five objects with follow-up near-infrared spectroscopy (Wu et al., 2018), roughly half of HotDOGs are consistent with a type 2 optical classification – i.e., with permitted lines of the same shape as the forbidden lines and no evidence for a separate broad-line region (Reyes et al., 2008). In contrast, the majority of ERQs show evidence for a broad-line region in their optical spectra (Zakamska et al., 2016; Perrotta et al., 2019). As a result, Fan et al. (2016) and Farrah et al. (2017) assume that the nucleus does not contribute at rest-frame optical wavelengths and do not perform PSF subtraction.

Glikman et al. (2015) analyze HST images of radio+infrared selected red quasars with bolometric luminosities $L_{\rm bol}=10^{47.3-48.3}$ erg/sec at $z=1.7-2.3$. Their observations are obtained with F160W and a shorter wavelength filter (F105W or F125W) over one or two orbits. They conduct PSF+Sersic decomposition, as well as merger classification and asymmetry measurements, and our observation and analysis methods are most directly comparable to this sample.

Finally, Hilbert et al. (2016) present the initial results of their study of hosts of radio-loud type 2 and type 1 active nuclei with F140W at $z\ga 1$. The redshifts of their objects are lower than ours, but they have three radio galaxies at $z=1.8-2.5$ which we include in our qualitative comparisons: 3C257, 3C326.1 and 3C454.1. In these three cases we assume that the nucleus does not contribute to the integrated photometry due to the type 2 spectrum in the rest-frame optical (e.g., Nesvadba et al. 2017 for 3C257).

4.5 Comparison of our sample with others

Despite the uncertainties in the fitting parameters, our $B$-band luminosities seem robust because the results of the parametric (Sersic) fits agree well with the results of the non-parametric (aperture photometry) estimates. Given that the stellar mass-to-light ratios are not known either for our sample or for any of the comparison ones, we restrict ourselves to direct comparison of host luminosities, not stellar mass. Furthermore, we calculate k-corrections as described in Section 3.4 for all comparison datasets, so that we can compare $L_{B}$ among the different samples.

Three of the comparison samples have $B$-band host galaxy luminosity distributions that are statistically consistent with those of ERQs in the sense of the Kolmogorov-Smirnov test. Aperture fluxes of HotDOGs measured by Farrah et al. (2017) correspond to $\log(L_{B}/L_{\odot})=11.1,11.2\pm 0.5$ (median, mean and sample standard deviation). These luminosities may be upper limits, since in the absence of near-infrared spectroscopy it is not yet known whether the quasars make an appreciable contribution to this luminosity and it is not clear whether the apertures exclude close companions. For their sample of type 1 quasars, Mechtley et al. (2016) present $V$-band luminosities which we k-correct to $B$-band by $-0.19$ dex. For their sample we find host luminosities $\log(L_{B}/L_{\odot})=11.1,11.1\pm 0.3$, again quite consistent with those of ERQs. The three radio galaxies in Hilbert et al. (2016) that have redshifts from $z=2.5$ to $z=1.8$ have luminosities $\log(L_{B}/L_{\odot})=11.6,11.3\pm 0.47$, somewhat larger than, but statistically consistent with, those of ERQs given the low number counts. This latter comparison comes with a caveat that 1 Gyr elapsed between these two epochs, so the two populations may be quite distinct, especially if the hosts of our targets average star formation rates of 100 $M_{\odot}$/year or more.

We estimate that the median stellar mass of the host galaxies of our targets ranges between $10^{10.4}-10^{11.4}M_{\odot}$ by assuming a range of plausible mass-to-light ratios (Maraston, 2005). The stellar masses could be higher if the host galaxies are extremely obscured on several kpc scales. Since all the comparison studies listed above derive their estimated host stellar masses from rest-frame $B$ or $V$-band data similar to ours, their estimates are subject to similar types of systematic uncertainties. Interestingly, our range of stellar masses is consistent with the stellar masses of high-redshift radio galaxies derived using completely different methods (Seymour et al., 2007).

The hosts of radio-selected red quasars from Glikman et al. (2015) are marginally brighter than those of ERQs, with $\log(L_{B}/L_{\odot})=11.6,11.5\pm 0.6$. However, out of the 10 radio-selected red quasars in the Glikman et al. (2015) sample, only five hosts are well-measured using PSF Sersic fits, and therefore the average $L_{B}$ quoted above is likely overestimated. The surprising property of this sample is that in about half of the sources, the quasar host galaxies are fainter than the clearly detected companions, whereas in our sample there is not a single object where the companion is brighter than the host.

We do not find strong evidence that the host galaxies of ERQs or type 2 quasar candidates are associated with ongoing major mergers. The fraction of objects with bright nearby companions in either one of our two subsamples is only 20$-$30%, short of theoretical expectations from numerical simulations (Steinborn et al., 2018; McAlpine et al., 2018), especially for the most luminous quasars at this redshift (Hickox et al., 2014). In observations, some recent works (Glikman et al., 2015; Fan et al., 2016; Hilbert et al., 2016) find a 80$-$100% merger incidence in luminous obscured populations. For HotDOGs, Farrah et al. (2017) argue that the high merger incidence derived by Fan et al. (2016) may be over-estimated. Visually, about 40$-$50% of Farrah et al. (2017) sources would likely be classified by us as ongoing major mergers, which (given the small number statistics) may be consistent with our sample. Their asymmetries are even lower than ours, though they are not directly comparable due to PSF treatment differences. In contrast, high-redshift radio-loud samples may show genuinely higher merger fraction than our sources, as we discuss in Section 4.6.

Submm observations of a HotDOG by Diaz-Santos et al. (2018) and of a luminous type 1 quasar by Bischetti et al. (2018) reveal a multitude of gas companions and tidal tails in host galaxies which are not detectable in the rest-frame optical data (similar techniques are also used in uncovering the environments of extremely high-redshift quasars, Decarli et al. 2017). Perhaps in ERQs merger activity can similarly be hidden in rest-frame optical observations. Alternatively our targets could be associated with a late stage of the merger. They are less obscured than HotDOGs, which may or may not be an orientation effect, so the marginally lower merger incidence rate in our sample compared to HotDOGs may indicate that ERQs are more likely to occur in a later phase (Hamann et al., 2017). Finally, perhaps extremely luminous quasar activity does not have to be associated with mergers at all and might result from processes within the galaxy (Shlosman et al., 1989). Whether relevant amounts of gas can be supplied by such processes to power the likely super-Eddington ERQs should be explored with numerical simulations. Type 2 quasar candidates in our sample do not show luminosities that are high enough to suggest super-Eddington accretion, and they may be suitable candidates for such galactic processes.

4.6 Comparison with other samples: red radio quasars and 3C sources

The merger incidence of ERQs and type 2 quasar candidates is not high and is not statistically distinguishable from that of Farrah et al. (2017) or Mechtley et al. (2016) with the current data, and Mechtley et al. (2016) in turn suggest that the incidence of mergers in their quasar sample is indistinguishable from that of high-redshift field galaxies of the same mass. Either the types of quasar activity we target in this paper are not predominantly associated with ongoing major mergers, or they may be associated with a post-merger phase we cannot identify using our data. The sample by Glikman et al. (2015) shows qualitatively different behavior. These objects are frequently found in ongoing mergers of highly disturbed hosts and are often hosted by the less massive companion in a pair, and their asymmetries (measured exactly using the same method as ours) are in the range $0.5-0.8$, marginally higher than our values, with host galaxy images visually more disturbed than ours. However, just like Glikman et al. (2015), we do not regard our asymmetries as reliably measured and are including them in our discussion mostly to emphasize the difficulties of quantitative morphology measurements in the presence of a bright central source.

These targets are similar to our ERQs in their bolometric luminosities and red colours, but have somewhat lower obscuration (Glikman et al., 2017). So why are we reaching different conclusions regarding the morphologies of our ERQ sample? An interesting possibility is that the host morphology may be tied to the radio properties. The Glikman et al. (2015) sources are radio-selected, and their median radio luminosity measured from FIRST (Becker et al., 1995) survey fluxes and k-corrected to the rest-frame 1.4 GHz using $\alpha=-0.7$ is $\nu L_{\nu}=10^{41.9}$ erg/sec. In contrast, neither ERQs nor type 2 quasar candidates are radio-selected. The median luminosity of the ten ERQs presented here is $\nu L_{\nu}$[1.4GHz]$=10^{40.9}$ erg/sec (Hwang et al., 2018), an order of magnitude lower, and type 2 quasar candidates are typically even less radio luminous (Alexandroff et al., 2016).

There are other indications that radio loudness may be related to host properties in high-redshift populations, with very high merger fractions reported for radio-loud quasars by Ivison et al. (2012) and Noirot et al. (2018). Chiaberge et al. (2015) find that at $z\simeq 1$ obscured (type 2) 3C radio galaxies are predominantly associated with mergers, whereas radio-quiet type 2 quasars are not any more likely to be in mergers than are non-active galaxies. Hilbert et al. (2016) present a larger sample of obscured and unobscured radio-loud quasars at $1<z<2.5$. They find that all obscured radio-loud hosts, where PSF subtraction is not necessary for evaluating the morphology, are highly disturbed, have multiple companions and appear to be associated with mergers. Of the three objects at the same redshifts as our targets, two would be classified by us as ongoing major mergers and the third has several faint companions. The analysis of the unobscured radio-loud quasars from this sample (Chiaberge et al., 2018) which uses PSF subtraction finds a strong association with mergers among these objects as well. More indirectly, radio-loud sources, including those at redshifts comparable to those of our sample, are found in more dense large-scale environments (Wylezalek et al., 2013) than comparable radio-quiet objects, which may naturally result in a higher merger fraction.

We conclude that radio luminosity is an important parameter to control for in future studies of quasar hosts at high redshift. As argued by Hwang et al. (2018), the radio luminosity of ERQs may be on the upper end of what may be considered a “radio-quiet” population not associated with powerful jets. Therefore, even though the radio luminosity difference between ERQs and Glikman et al. (2015) is only an order of magnitude, we speculate that the two samples may have dramatically different mechanisms of producing radio emission, with jets – relativistic, collimated outflows – in the Glikman et al. (2015) sources and winds – non-relativistic, non-collimated outflows – in ERQs. The 3C objects from Hilbert et al. (2016) and (Chiaberge et al., 2018) are on the extreme bright end of the radio-loud distribution and are clearly jet-dominated. Correspondingly, it is possible that either the parameters of the supermassive black holes (e.g., spins) or the accretion parameters might be quite different in the populations with and without jets, despite having similar bolometric luminosities. Thus it may not be surprising that the host properties and the larger scale environments may be different as well.

Another possibility for the difference between our sample and those of Glikman et al. (2015) and Hilbert et al. (2016) is that the merger fraction may be higher in populations selected in the infrared and/or in the radio without regard to optical (rest-frame UV) fluxes. Such selections may naturally result in a bias toward populations with higher galaxy extinction and therefore toward a higher merger fraction. In support of this possibility is the conclusion reached by Glikman et al. (2017) that the obscuration in their sample arises on the scales of the host galaxy, whereas in ERQs it is likely circumnuclear (Goulding et al., 2018b). However, this hypothesis fails to account for the modest merger fraction of HotDOGs (selected based on the mid-infrared alone) which is comparable to that seen in our sample (Farrah et al., 2017). Therefore, the connection between radio loudness and host merger activity at high redshift is well worth exploring with future observations with HST and JWST.

5 Conclusions

In this paper we present HST observations of ten ERQs and six type 2 quasar candidates at $2<z<3$ – the peak epoch of quasar activity and galaxy formation. ERQs are extremely luminous (up to $L_{\rm bol}\la 10^{48}$ erg/sec) near- or super-Eddington accretors, with unprecedented ionized gas outflow activity seen both on small and large scales. No analogs to these objects have been found in the low-redshift universe. ERQs are ideal candidates for being observed during the early short-lived feedback phase of quasars when they are still dust enshrouded but driving extremely powerful outflows that expel the gas and dust, disrupt star formation, and halt mass assembly in the host galaxies. While this evolution is expected to be triggered by galaxy mergers and interactions (Hopkins et al., 2006), there is not yet a definitive observation of this connection during the epoch of peak galaxy formation and quasar activity. Type 2 quasar candidates are less powerful, show moderate amount of ionized wind activity in their [OIII] emission lines, and resemble low-redshift type 2 quasars, except perhaps at lower levels of obscuration due to their original selection at rest-frame ultraviolet wavelengths.

In our HST programmes, we image both samples in the ACS F814W and WFC3 F160W (or F140W) bands. We use a careful PSF subtraction procedure on the WFC3 data to reveal the host galaxy emission in almost every object. We also perform joint PSF+Sersic fits to determine the structural parameters of the host galaxies.

The median rest-frame $B$-band luminosities of the hosts are $10^{11.2}L_{\odot}$ and $10^{10.6}L_{\odot}$ for ERQs and for type 2 quasar candidates, respectively, with the latter value similar to the $L^{*}$ luminosity at this redshift. Despite uncertainties in PSF subtraction and difficulty of joint PSF+Sersic fits, there is good agreement between several different measurements of host luminosities, and we estimate the uncertainties in luminosity of each object at $\la 0.25$ dex. These luminosities are similar to the $B$-band luminosities found for other luminous quasar hosts using similar data. The structural parameters, such as effective radii and Sersic indices, are more uncertain and more sensitive to the exact fitting procedure. Most fits indicate that the galaxies in our sample are early-type, with Sersic indices between 2 and 3 and a median effective radius just under 5 kpc; however, Sersic indices have especially large estimated uncertainties ($\sigma(n_{s})\simeq 1$) and are not yet well constrained. These properties are similar to those of massive high-redshift galaxies in van Dokkum et al. (2010), except our hosts are somewhat less compact – with $R_{50}$ of 5 kpc instead of $<$3 kpc – but this may well be within the large uncertainty on $R_{50}$.

Only $20-30$% of ERQs and type 2 quasar candidates are in ongoing major mergers with a flux ratio greater than 1:4 and with components separated by $<$2.5″. This is in contrast to the $80-100$% merger fraction expected for quasars at this luminosity and redshift and found among the high-redshift radio-loud population. We have discussed three possibilities for this unexpected result. (i) Merger signatures can be hidden at rest-frame optical wavelengths in high-redshift dusty galaxies and may be revealed by future submm observations. (ii) Objects in our sample might not be associated with mergers at all, and may instead be powered by gas instabilities within galactic discs. This possibility is more likely to apply to the type 2 quasar candidates which are less luminous than ERQs and are thus less likely to be super-Eddington accretors. Numerical simulations are necessary to determine whether super-Eddington accretion of ERQs can result from galactic processes in the absence of mergers. (iii) Objects in our sample might be associated with late-stage mergers that cannot be identified using our HST data due to the faintness of the tidal features. Higher sensitivity JWST observations can be used to probe this scenario.

Acknowledgements

Based on observations made with the NASA/ESA HST, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programmes GO-14608 and GO-13014. Support for these programmes was provided by NASA through grants HST-GO-14608 and HST-GO-13014 from the STScI.

NLZ acknowledges support by the Catalyst Award at Johns Hopkins University and by the Deborah Lunder and Alan Ezekowitz Founders’ Circle Membership at the Institute for Advanced Study where part of this work was conducted. WNB acknowledges support by the NASA ADP program and NSF grant AST-1516784. GL is supported by the National Thousand Young Talents Program of China and acknowledges the grant from the National Natural Science Foundation of China (No. 11673020 and No. 11421303) and the Ministry of Science and Technology of China (National Key Program for Science and Technology Research and Development, No. 2016YFA0400700). NPR acknowledges support from the STFC and the Ernest Rutherford Fellowship scheme.