[OIII] Emission Line Properties in a New Sample of Heavily Reddened Quasars at $z>2$

TL;DR

This study analyzes the [OIII] emission line properties in a new sample of heavily reddened quasars at redshifts greater than 2, revealing strong outflows and suggesting dust is located on galactic scales.

Contribution

First comparison of [OIII] line properties in high-luminosity obscured quasars to unobscured ones at similar redshifts and luminosities, providing insights into outflows and dust location.

Findings

Broad wings in [OIII] lines indicate strong outflows.

No significant kinematic difference between reddened and unobscured quasars.

Dust likely resides on galactic scales outside the [OIII] emitting region.

Abstract

We present VLT-SINFONI near infra-red spectra of 26 new heavily reddened quasar candidates selected from the UKIDSS-LAS, VISTA VHS and VIKING imaging surveys. This new sample extends our reddened quasar search to both brighter and fainter luminosities. 25 of the 26 candidates are confirmed to be broad line quasars with redshifts $0.7<z<2.6$ and dust extinctions $0.5<E(B-V)<3.0$. Combining with previously identified samples, we study the H$\alpha$, H$\beta$ and [OIII] emission line properties in 22 heavily reddened quasars with $L_{bol}\approx 10^{47}$ erg/s and $z>2$. We present the first comparison of the [OIII] line properties in high luminosity obscured quasars to a large sample of 111 unobscured quasars in the same luminosity and redshift range. Broad wings extending to velocities of 2500 km/s are seen in the [OIII] emission line profiles of our reddened quasars, suggesting that strong outflows are affecting the ionised gas kinematics. However, we find no significant difference between the kinematics of the [OIII] emission in reddened and unobscured quasars when the two samples are matched in luminosity and redshift. Our results are consistent with a model where quasar-driven outflows persist for some time after the obscuring dust has been cleared from along the line of sight. Assuming the amount of ionised gas in obscured and unobscured quasars is similar, we use the equivalent width distribution of the [OIII] emission to constrain the location of the obscuring dust in our reddened quasars. We find that the dust is most likely to be located on galactic scales, outside the [OIII] emitting region.

Tables

Table 1. Table 1: Observation log for our 26 targets. Exposure times are given in seconds and the full width at half maximum of the point spread function (i.e., the size of the seeing disc) in arcseconds. Redshifts are taken from the broad component of the fit to H  α 𝛼 \alpha , where present, otherwise from the Paschen lines. J ​ H ​ K 𝐽 𝐻 𝐾 JHK photometry are taken from UKIDSS-LAS for ULAS sources, and J ​ H ​ K s 𝐽 𝐻 subscript 𝐾 𝑠 JHK_{s} photometry from VISTA for VHS and VIK sources. Spectra are shown in Fig.  8 .

Name	RA	Dec.	$J_{\mathrm{Vega}}$	$H_{\mathrm{Vega}}$	$K_{\mathrm{Vega}}$	Redshift	Observation date	Exp. time (s)	Seeing (″)
VHS J1117-1528	11:17:03.74	-15:28:29.6	20.75	18.44	16.94	2.428	2015-04-06	2400	1.20
VHS J1122-1919	11:22:24.43	-19:19:17.4	19.58	-	15.96	2.464	2015-04-06	600	1.52
ULAS J1216-0313	12:16:31.78	-03:13:35.0	-	18.17	16.20	2.574	2015-04-06	600	1.38
VHS J1301-1624	13:01:31.32	-16:24:54.0	20.12	-	16.39	2.138	2015-04-07	1200	1.06
ULAS J1415+0836	14:15:26.74	+08:36:15.8	-	18.67	16.64	1.120	2015-06-05	1200	1.88
VIK J2205-3132	22:05:13.68	-31:32:02.4	19.93	18.46	17.33	2.307	2015-06-15	2400	0.95
VIK J2214-3100	22:14:05.40	-31:00:38.2	20.25	18.67	16.84	1.069	2015-06-22	2400	0.75
VHS J2222-3915	22:22:10.65	-39:15:49.3	18.77	-	14.60	-	2015-04-04	400	0.85
VIK J2228-3205	22:28:44.98	-32:05:38.8	21.00	19.69	18.40	2.364	2015-07-07	1200	1.22
							2015-07-02	2400	0.96
VIK J2230-2956	22:30:09.36	-29:56:21.5	20.93	18.98	17.47	1.319	2015-07-07	1200	1.36
							2015-07-02	2400	0.87
VIK J2232-2844	22:32:36.96	-28:44:39.5	20.02	18.57	16.93	2.292	2015-08-07	2400	0.99
VIK J2238-2836	22:38:23.95	-28:36:40.7	19.56	18.14	17.03	1.231	2015-07-28	2400	1.15
VIK J2241-3006	22:41:12.41	-30:06:42.5	19.56	18.27	16.84	0.720	2015-07-21	2400	0.89
VIK J2243-3504	22:43:48.96	-35:04:39.4	19.78	18.01	16.41	2.085	2015-07-20	1200	1.12
VIK J2245-3516	22:45:58.80	-35:16:58.4	20.30	18.81	17.17	1.335	2015-08-07	2400	1.31
							2015-07-02	600	1.19
VIK J2251-3433	22:51:17.90	-34:33:50.0	21.10	18.86	18.04	1.693	2015-08-07	2400	1.31
							2015-08-11	1200	1.64
VIK J2256-3114	22:56:07.97	-31:14:25.8	20.01	18.55	17.25	2.329	2015-07-21	2400	0.80
VIK J2258-3219	22:58:30.10	-32:19:11.6	18.96	17.69	16.45	0.879	2015-05-10	1200	1.67
VIK J2306-3050	23:06:50.23	-30:50:34.1	20.72	19.25	18.07	1.060	2015-07-07	1200	1.86
							2015-07-07	2400	1.31
VIK J2309-3433	23:09:49.90	-34:33:43.6	20.77	19.50	18.26	2.159	2015-07-19	1200	0.74
							2015-07-19	2400	0.74
ULAS J2312+0454	23:12:31.99	+04:54:20.9	20.07	18.67	16.72	0.700	2015-08-01	1200	1.20
VIK J2313-2904	23:13:51.96	-29:04:52.7	20.22	18.00	17.60	1.851	2015-07-10	1200	0.85
							2015-07-02	2400	1.10
VIK J2314-3459	23:14:16.87	-34:59:47.0	19.89	18.50	17.36	2.325	2015-07-20	2400	0.94
VIK J2323-3222	23:23:50.45	-32:22:16.7	20.12	18.44	16.83	2.191	2015-05-10	1200	1.13
VIK J2350-3019	23:50:45.50	-30:19:53.4	20.12	18.81	17.58	2.324	2015-07-19	2400	0.65
							2015-06-02	1200	1.24
VIK J2357-3024	23:57:45.84	-30:24:30.6	20.33	18.75	17.58	1.129	2015-07-02	2400	0.94
							2015-07-30	1200	0.97

Table 2. Table 2: Derived properties from our SINFONI observations of the 12 new quasars at z < 2 𝑧 2 z<2 . Black hole mass estimates have uncertainties of ∼ 0.4 similar-to absent 0.4 \sim 0.4 dex.

Name	Redshift	${\left(H-K\right)}_{\mathrm{AB},\mathrm{SINFONI}}$	$E\left(B-V\right)$	${\text{Log}}_{10}\left({\text{M}}_{BH}/M_{\odot }\right)$
ULAS J1415+0836	1.120	1.73	3.0	9.03
VIK J2214-3100	1.069	1.48	2.5	8.80
VIK J2230-2956	1.319	1.04	2.0	8.63
VIK J2238-2836	1.231	0.87	1.4	8.99
VIK J2241-3006	0.720	1.20	2.0	7.56
VIK J2245-3516	1.335	1.66	3.0	9.59
VIK J2251-3433	1.693	1.05	1.9	9.23
VIK J2258-3219	0.879	0.84	1.3	7.52
VIK J2306-3050	1.060	1.07	1.8	8.34
ULAS J2312+0454	0.700	1.36	2.4	8.24
VIK J2313-2904	1.851	1.10	1.4	9.58
VIK J2357-3024	1.129	1.36	2.2	8.66

Table 3. Table 3: Derived properties for the 22 quasars at z > 2 𝑧 2 z>2 for which we have coverage of H  α 𝛼 \alpha , H  β 𝛽 \beta , and [ O  iii ]. Redshifts are taken from the broad component of the fit to the H  α 𝛼 \alpha line. The E ​ ( B − V ) 𝐸 𝐵 𝑉 E(B-V) is estimated from the near infra-red photometry, and has an uncertainty of ∼ 0.15 similar-to absent 0.15 \sim 0.15 magnitudes. Rest frame 5100 Å luminosities have been corrected for dust extinction using the E ​ ( B − V ) 𝐸 𝐵 𝑉 E(B-V) , giving rise to an uncertainty of ∼ 0.2 similar-to absent 0.2 \sim 0.2 dex. Black hole masses are derived using the prescription from Vestergaard & Peterson ( 2006 ) , and have an uncertainty of ∼ 0.4 similar-to absent 0.4 \sim 0.4 dex. Eddington ratios are estimated assuming L = 8 × L 5100 𝐿 8 subscript 𝐿 5100 L=8\times L_{5100} , and have an associated uncertainty of a factor of 4. H  α 𝛼 \alpha equivalent widths are given in the rest frame. [ O  iii ]  λ 𝜆 \lambda 5008 line luminosities are given in the observed frame with no correction for extinction. 80 per cent velocity widths have been corrected for instrumental broadening of σ = 200 ​ km ​ s − 1 𝜎 200 km superscript s 1 \sigma=200\hbox{$\rm\thinspace km\rm\thinspace s^{-1}$} . *Indicates that object is consistent with having no [ O  iii ] emission, and an upper bound on the line luminosity has been estimated from a template fit.

Name	Redshift	$E\left(B-V\right)$	${\text{Log}}_{10}\left({\text{M}}_{BH}\right)$	${\text{Log}}_{10}\left(L_{5100}\right)$	$L/L_{\mathrm{Edd}}$	EW(H $\alpha$)	${\text{Log}}_{10}\left(L_{\left[\text{O }\text{iii}\right]}\right)$	[O iii] $w_{80}$
		[magnitudes]	[$M_{\odot }$]	[$\mathrm{erg}{s}^{-1}$]		[Å]	[$\mathrm{erg}{s}^{-1}$]	[$\mathrm{km}{s}^{-1}$]
ULAS J0016-0038	2.176	0.55	9.20	45.95	0.36	883	43.22$\pm$0.07	730$\pm$131
ULAS J0041-0021	2.510	0.81	9.52	46.52	0.64	281	43.18$\pm$0.06	868$\pm$101
ULAS J0141+0101	2.562	0.58	9.43	45.97	0.22	490	43.79$\pm$0.05	2091$\pm$83
ULAS J0144-0014	2.504	0.68	9.58	46.40	0.42	286	-	-
ULAS J0144+0036	2.281	0.74	9.65	46.31	0.29	377	43.31$\pm$0.06	1880$\pm$89
ULAS J0221-0019	2.248	0.58	8.86	45.92	0.73	656	43.24$\pm$0.06	1093$\pm$108
VHS J1117-1528	2.428	1.25	9.52	46.65	0.87	342	43.16$\pm$0.06	2277$\pm$227
VHS J1122-1919	2.464	1.09	9.24	46.82	2.39	513	43.81$\pm$0.05	1775$\pm$130
ULAS J1216-0313	2.574	1.62	10.07	47.51	1.73	229	*	-
ULAS J1234+0907	2.502	1.94	10.48	47.64	0.91	432	*	-
VHS J1301-1624	2.138	1.24	9.40	46.59	0.99	364	43.41$\pm$0.06	2898$\pm$241
ULAS J2200+0056	2.541	0.55	9.10	46.47	1.49	461	44.15$\pm$0.05	1729$\pm$76
VIK J2205-3132	2.307	0.58	8.79	45.60	0.41	485	42.71$\pm$0.07	772$\pm$96
ULAS J2224-0015	2.219	0.61	8.98	46.04	0.73	624	43.20$\pm$0.07	1292$\pm$233
VIK J2228-3205	2.364	0.57	8.48	45.19	0.32	560	43.24$\pm$0.05	624$\pm$28
VIK J2232-2844	2.292	0.65	9.02	45.43	0.16	762	42.47$\pm$0.09	3161$\pm$1062
VIK J2243-3504	2.085	1.05	8.99	46.27	1.21	445	43.68$\pm$0.05	1710$\pm$75
VIK J2256-3114	2.329	0.69	9.16	45.77	0.26	394	43.18$\pm$0.05	1702$\pm$103
VIK J2309-3433	2.159	0.52	9.39	45.03	0.03	647	42.54$\pm$0.06	315$\pm$51
VIK J2314-3459	2.325	0.53	8.41	45.21	0.40	414	43.32$\pm$0.05	2572$\pm$341
VIK J2323-3222	2.191	0.99	9.09	46.18	0.79	475	42.98$\pm$0.06	979$\pm$273
VIK J2350-3019	2.324	0.58	8.58	45.57	0.62	337	43.03$\pm$0.05	1059$\pm$59

Equations

$$A_{\lambda }=k(\lambda )\times E(B-V)$$

$${\mathrm{FWHM}}_{\mathrm{H}\beta }=1.23\times 10^3{\left(\frac{{\mathrm{FWHM}}_{\mathrm{H}\alpha }}{10^3\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{0.97}$$

$${\mathrm{M}}_{\mathrm{BH}}=10^{6.91}{\left(\frac{{\mathrm{L}}_{5100}}{10^{44}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{0.5}{\left(\frac{{\mathrm{FWHM}}_{\mathrm{H}\beta }}{10^3\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2{\mathrm{M}}_{\odot }$$

$${\mathrm{M}}_{\mathrm{BH}}=2\times 10^6{\left(\frac{{\mathrm{L}}_{\mathrm{H}\alpha }}{10^{42}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{0.55}{\left(\frac{{\mathrm{FWHM}}_{\mathrm{H}\alpha }}{10^3\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{2.06}{\mathrm{M}}_{\odot }$$

$${\mathrm{M}}_{\mathrm{BH}}=10^{7.13}{\left(\frac{{\mathrm{L}}_{\mathrm{Pa}\beta }}{10^{42}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{0.48}{\left(\frac{{\mathrm{FWHM}}_{\mathrm{Pa}\beta }}{10^3\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2{\mathrm{M}}_{\odot }$$

Full text

††pubyear: 2019††pagerange: [O iii] Emission Line Properties in a New Sample of Heavily Reddened Quasars at $z>2$–B

1 Introduction

Feedback from active galactic nuclei (AGN) is used in the study of galaxy evolution to explain many observables, such as the shape of the galaxy luminosity function (e.g. Vogelsberger et al., 2013; Sijacki et al., 2015) and correlations between the mass of the central supermassive black hole (SMBH) and host galaxy properties (Magorrian et al., 1998; McConnell & Ma, 2013). Evidence for AGN feedback has been found through observations of different gas phases and at different wavelengths (e.g. Fabian, 2012; Harrison, 2017).

The most massive galaxies in today’s universe must have undergone at least one major merger, and possibly many more, in order to build up their present day masses; thus major mergers are an important component in many models of galaxy evolution. Such merger-driven models normally predict an reddened quasar phase: as the galaxy merger triggers a burst of intense star formation (which in turn produces large amounts of reddening dust), at the same time gas is dynamically shocked and falls onto the central SMBH, triggering highly luminous AGN activity (Sanders et al., 1988; Di Matteo et al., 2005; Hopkins et al., 2008; Narayanan et al., 2010).

This proposed model for quasar activity and galaxy evolution is reliant on a ‘blow-out’ phase, where powerful small scale outflows from the inner regions of the active nucleus couple to gas and dust, transferring energy and momentum and driving large scale outflows which then disrupt the gas in the host galaxy, shutting down star formation and clearing obscuring material away from the line of sight to reveal a luminous blue quasar. Such outflows are likely to be driven by some combination of gas pressure and radiation pressure acting on dust, neutral hydrogen and atomic line transitions, and disk-driven magnetocentrifugal winds (e.g. Konigl & Kartje, 1994; Murray et al., 1995; Ishibashi & Fabian, 2015; Ishibashi et al., 2017), although the dominant driving mechanism is still a topic of debate. If this model is correct, one might therefore expect to see evidence for more powerful outflows in observations of luminous reddened quasars than in unobscured quasars, while the obscuring material is in the process of being cleared from the line of sight.

One probe of large scale outflows in quasars is the forbidden [O iii] $\lambda\lambda$4960,5008 doublet transition, which traces ionised gas in regions of low electron density ($n_{e}\sim 10^{3}\rm\thinspace cm^{-3}$; Baskin & Laor, 2005). [O iii] is present as a strong emission line in many active galaxies, and the shape of the line can be used as a robust estimator of the kinematics of the emitting ionised gas. Spatially resolved observations of local (i.e. $z<0.5$, and therefore less luminous) active galaxies have found ionised gas moving with near-escape velocities on scales of $0.5-15\rm\thinspace kpc$ in most objects (Greene et al., 2011; Greene et al., 2012; Liu et al., 2013; Harrison et al., 2014; Tadhunter et al., 2018), although different observers use different definitions of spatial extent depending on whether they use spatially resolved kinematic data, narrow-band imaging or other methods (Karouzos et al., 2016; Baron & Netzer, 2019).

While the classic narrow [O iii] emission line profile can be used to define the systemic redshift, more luminous AGN often have a broader (v\sim 1000\hbox{\rm\thinspace km\rm\thinspace s^{-1}}) blueshifted component to the [O iii] line, which is too broad to be tracing gas in dynamical equilibrium with the host galaxy. If present, broad blue [O iii] emission is therefore usually interpreted as evidence for quasar-driven outflows on kiloparsec scales, driving the emitting gas along the line of sight towards the observer. If these outflows are indeed occurring on such large scales, then they would be capable of transferring significant amounts of kinetic energy from the active nucleus back into the host galaxy (Liu et al., 2013; Harrison et al., 2018; Coatman et al., 2019).

Using the [O iii] emission line, several authors have recently investigated the ionised gas kinematics in AGN across a range of redshifts, luminosities, and levels of obscuration. DiPompeo et al. (2018) find that at low redshifts ($z<0.4;L_{\rm bol}\approx 10^{46}\rm\thinspace erg\rm\thinspace s^{-1}$), AGN with redder $u-W3$ colours have broader, more blueshifted [O iii] on average. This difference is not seen between spectroscopically classified type 1 and type 2 quasars in the same sample, suggesting that the difference in ionised gas kinematics is not purely due to orientation, but that evolutionary effects might also be at work.

Harrison et al. (2016) find that the distribution of ionised gas velocities in an X-ray selected (L_{2-10\,\rm{keV}}>10^{42}\hbox{\rm\thinspace erg\rm\thinspace s^{-1},}) sample of AGN at $1.1<z<1.7$ is consistent with that of a luminosity-matched sample of $z<0.4$ AGN, suggesting that the most powerful ionised gas outflows are driven by AGN activity at all redshifts. Brusa et al. (2015) find that quasars selected to have bright X-ray fluxes and red optical–to–infra-red colours at $z\approx 1.5$ and L_{\rm bol}\approx 10^{45.5}\hbox{\rm\thinspace erg\rm\thinspace s^{-1},} have [O iii] lines which are broader than Sloan Digital Sky Survey (SDSS) type 2 AGN at similar [O iii] luminosities.

At higher redshifts ($z\approx 2.5$; $L_{\rm bol}\approx 10^{47}\rm\thinspace erg\rm\thinspace s^{-1}$), Zakamska et al. (2016) find broad, strongly blueshifted [O iii] in four of the most extreme quasars in the SDSS. These objects were selected to have strong C iv emission and red optical–to–infra-red colours. Zakamska et al. argue that these objects may represent a ‘blow-out’ phase, as the emitting gas cannot be confined to the gravitational potential of the galaxy. Such ‘Extremely Red Quasars’ have extinctions $E(B-V)$ $\approx 0.3$ (Hamann et al., 2017).

The evidence so far therefore suggests that red and/or reddened quasars generally have stronger ionised gas outflows, and has often been interpreted as showing that the red selection (i.e. the selection of objects with more obscuring dust) is responsible for the observed outflow properties. However, there has been no systematic, luminosity-matched comparison of reddened and unobscured quasars at $z>2$, corresponding to the epoch of peak star formation and peak AGN activity (Madau & Dickinson, 2014; Kulkarni et al., 2018).

Using near infra-red data, Banerji et al. (2012); Banerji et al. (2013, 2015) have discovered a population of heavily dust-reddened ($E(B-V)\approx 1$), broad line (i.e., spectroscopic type 1) quasars at redshifts $z\approx 2$. These quasars have been selected to have extinctions similar to those in starburst galaxies, as one would expect following an intense burst of star formation triggered by a major merger. At the same time they show broad emission lines, proving that the obscuration is not solely due to simple orientation effects and also allowing their black hole masses to be estimated from virial methods. Four of these quasars have been followed up with the Atacama Large Millimeter/sub-millimeter Array (ALMA; Banerji et al., 2017; Banerji et al., 2018), finding high levels of both cold dust and molecular gas in their host galaxies, while Wethers et al. (2018) exploit the obscuration of the quasar light in the rest frame UV to investigate the host galaxy emission, finding evidence for high rates of ongoing star formation.

In this work, we present the first direct comparison of the properties of the ionised gas emission in the unobscured and heavily dust-reddened quasar populations at L_{\rm bol}\approx 10^{47}\hbox{\rm\thinspace erg\rm\thinspace s^{-1},} and $z>2$. We make use of the catalogue of near infra-red spectroscopic observations from Coatman et al. (2019), which describes the ionised gas kinematics in a large ($>100$) sample of unobscured quasars in this luminosity and redshift range for the first time. Quantifying the properties of the ionised gas emission in reddened and unobscured quasars allows us to compare the feedback mechanisms in the two populations, and simultaneously place constraints on the location of the obscuring material in the reddened quasars.

The structure of this paper is as follows. In Section 2 we present the photometric selection and spectroscopic follow-up of 25 new heavily reddened broad line quasars. By combining with previous studies, we then construct a sample of 22 heavily reddened quasars at redshift $z>2$ with spectral coverage of the [O iii] $\lambda$5008 emission line, and describe our fitting procedure to derive H $\alpha$ and [O iii] emission line properties in Section 3. We compare the ionised gas kinematics in our heavily reddened sample with a sample of unobscured quasars in the same redshift and luminosity range in Section 4 and discuss our results in Section 5.

We assume a flat $\Lambda$CDM cosmology throughout this work, with $\Omega_{m}=0.27$, $\Omega_{\Lambda}=0.73$, and \textrm{H}_{0}=71\hbox{\hbox{$\rm\thinspace km\rm\thinspace s^{-1}$}\rm\thinspace Mpc^{-1}}. All emission lines are identified with their wavelengths in vacuum in units of Ångströms. We avoid the inherent uncertainty involved in estimating bolometric luminosities by using the monochromatic continuum luminosity $\lambda L_{\lambda}$ at $\lambda=5100$ Å (hereafter $L_{5100}$) as a measure of the strength of the ionising radiation. When estimating Eddington ratios, we assume a bolometric correction $L_{\rm bol}=8\times L_{5100}$ (e.g. Nemmen & Brotherton, 2010; Runnoe et al., 2012). We correct for dust extinction using the quasar extinction law described in Wethers et al. (2018), which gives wavelength-dependent attenuations

[TABLE]

with $k(\rm{H}\,\beta)=3.57$, $k(5008\,$Å$)=3.45$, $k(5100\,$Å$)=3.38$, $k(V)=3.10$, $k(\rm{H}\,\alpha)=2.53$, $k(\rm{Pa}\,\gamma)=1.41$, and $k(\rm{Pa}\,\beta)=1.17$.

2 Data

2.1 Photometric sample selection

Our previous searches for heavily reddened quasars (Banerji et al., 2012; Banerji et al., 2013, 2015) have made use of near infra-red photometric observations in the $J$, $H$, and $K$-bands from the UKIDSS Large Area Survey (UKIDSS-LAS) and the VISTA Hemisphere Survey (VHS). In these previous studies, the wide-field near infra-red data was used to select extremely red quasar candidates e.g. with $(J-K)_{\rm{Vega}}>2.5$. Subsequent spectroscopic follow-up then confirmed that many of these candidates are quasars at $2<z<3$, with typical extinctions of $E(B-V)$ $\approx 0.8$ towards the quasar continuum. In Banerji et al. (2015), we constructed the first luminosity function for the reddened quasar population, which showed evidence that the reddened quasars outnumber unobscured quasars at the highest quasar luminosities, but that their number counts begin to decline as we approach more typical quasar luminosities around $L^{*}$. In this work, we extend the luminosity range covered by our reddened quasar sample by specifically searching for these quasars at the bright and faint ends of the luminosity function presented in Banerji et al. (2015).

2.1.1 Bright sample

For the selection of a very luminous sample of extremely red quasars, we make use of near infra-red imaging data over $\approx$6300 deg2 from the UKIDSS-LAS and VHS and mid infra-red imaging from WISE. The UKIDSS-LAS search is concentrated in three regions (i) 08h$<$RA$<$16h; Dec$<$20$\degr$ (1920 deg2) (ii) 20h$<$RA$<$00h; Dec$<$20$\degr$ (455 deg2) and (iii) 00h$<$RA$<$04h; Dec$<$20$\degr$ (755 deg2). The VHS search is similarly concentrated in three regions: (i) 10h$<$RA$<$16h; -35$<$Dec$<$0$\degr$ (1510 deg2) (ii) 20h$<$RA$<$00h; -25$<$Dec$<$-2°(670 deg2) and (iii) 20h$<$RA$<$00h; -65$<$Dec$<$-35°(980 deg2). All candidates were selected to be classified as point sources in the $K$-band and have $K_{\rm{Vega}}<17.0$. Where $H$-band data was available, we used a colour cut of $(H-K)_{\rm{Vega}}>1.9$ to select the reddest, most intrinsically luminous quasar candidates. This colour cut corresponds to $E(B-V)\gtrsim 1$ at $z=2$. For some regions in VHS, only $J$ and $K$-band data is available. We therefore used a colour cut of $(J-K)_{\rm{Vega}}>3.4$ in these regions, corresponding to the same extinction criterion, and also required that sources be fainter in the $Y$-band where available, i.e. $Y>J$. Finally, as in our previous searches for reddened quasars, we applied an additional colour cut based on the WISE bands of $(W1-W2)>0.85$. With these selection criteria, a total of 24 candidates were identified over the full 6300 deg2 area. Seven of these candidates were observed and presented in Banerji et al. (2012); Banerji et al. (2015). A further seven extremely red, bright quasar candidates are spectroscopically followed up and presented in this paper.

2.1.2 Faint sample

To extend the search to fainter intrinsic luminosities than is possible with the UKIDSS-LAS or VHS, we made use of deeper near infra-red data from the VISTA VIKING survey (Edge et al., 2013). VIKING is covering 1500 deg2 in southern and equatorial fields and targets were selected over a $\approx$210 deg2 region with 22h$<$RA$<$00h; -40$<$Dec$<$-25°. As in our previous work, we required all candidates to have $(J-K)_{\rm{Vega}}>2.5$ and $(W1-W2)>0.85$ but now going down to a fainter magnitude limit of $K_{\rm{Vega}}<18.4$. These initial colour cuts produced a large candidate list of 177 sources, 40 of which are classified as point sources in the $K$-band and are therefore likely to be dominated by the quasar light in the infra-red bands. As discussed extensively in Banerji et al. (2012) and Banerji et al. (2015), selecting point sources reduces contamination from low redshift, low luminosity quasars where the light in the $K$-band has significant contributions from the quasar host galaxy. To further reduce the number of candidates, we matched these 40 point sources to new optical data from the VST-KiDS survey and required all candidates to either be undetected in the KiDS $r$ and $i$-bands, or have very red colours of $(i-K)_{\rm AB}>3.5$, reducing our number of quasar candidates to 17. In order to fill our allocation of observing time, we also targeted two of the reddest extended $K$-band sources with no VST-KiDS detections and $(J-K)_{\rm Vega}>3$, producing a final sample of 19 candidates over the VIKING region.

2.2 New spectroscopic data

All 26 targets were followed up with ESO-SINFONI, a near infra-red Integral Field Unit installed at the Cassegrain focus of VLT-UT4. We use the SINFONI $H$+$K$ grating which provides a spectral resolution of $R\approx 1500$ in the wavelength range 1.45-2.45 $\mu\textrm{m}$. We use a field of view of 8″x 8″ which is nodded in four exposures with offsets of $\pm$1.5″ in RA and Dec, giving a usable field of view of 5″x 5″ in the reduced datacubes. For targets in the redshift range $2.0<z<2.6$, this setup yields coverage of the H $\alpha$ emission line in the $K$-band and both H $\beta$ and [O iii] in the $H$-band. For redshifts of $1.2<z<2.0$ we have coverage of H $\alpha$ in the $H$-band. Lower redshift objects can be identified from their near infra-red spectra using the hydrogen Paschen series.

The data were reduced using the easysinf package, as described in Williams et al. (2017). 1-dimensional spectra were then extracted from the datacubes using circular apertures with diameters in the range 0.6-1.3″ ($\approx 5$-10 kpc at $z=2$). The size of the aperture was chosen to maximise the signal-to-noise ratio (S/N) of the H $\alpha$ line in each spectrum. Choosing instead to maximise the S/N of the continuum or the [O iii] line does not significantly change the size of the aperture in any of the objects and does not affect the shape of the spectrum, meaning there is no evidence for spatially extended [O iii] emission on scales $\gtrsim 5\rm\thinspace kpc$ within the limited S/N of our observations. In the cases where a single object was observed across multiple observing blocks (OBs), spectra were extracted from each datacube before being combined, as different OBs were conducted under different seeing conditions, leading to different optimal aperture sizes. Noise arrays were extracted by considering the spaxel-spaxel variance in the field away from the source at each wavelength pixel. Absolute flux calibration was carried out for each spectrum using the $K$-band photometry from VISTA and UKIDSS.

Of the 26 targets which satisfied our selection criteria and were spectroscopically followed up, 25 are confirmed to be quasars in the redshift range $0.7<z<2.6$. One target (VHS J2222-3915) is found to have a red continuum with no emission features. These observations are summarised in Table 1 and spectra are shown in Fig. 8.

2.3 Existing spectroscopic data

13 of the 25 spectroscopically confirmed quasars from the new sample presented above are at $z>2$. We combine these 13 objects with nine quasars from Banerji et al. (2012) in the same redshift range and with the same wavelength coverage – three objects with VLT-SINFONI data and six objects with Gemini-GNIRS data – to obtain a total of 22 quasars with near infra-red spectra covering H $\alpha$, H $\beta$ and [O iii], which we refer to as the ‘high redshift sample’ and use when reporting results in the rest of this paper. We carry out our own fits to the emission lines in the previously published objects and hence derive parameters in a consistent manner for the whole sample.

3 Spectral fitting and black hole masses

3.1 High redshift sample

For the 22 quasars with $z>2$ in our high redshift sample, we have coverage of [O iii], which traces ionised gas in regions of low electron density. Measuring the shape and strength of the [O iii] emission therefore allows us to quantify the kinematics of the ionised gas in these objects. Using the kinematics of the Balmer emitting gas, we can also estimate black hole masses and Eddington ratios. We model the emission from H $\alpha$, H $\beta$ and [O iii] in order to derive more robust properties for these emission lines.

We now describe our modelling procedure for the H $\alpha$, H $\beta$ and [O iii] emission lines in our high redshift sample, which is very similar to that used by Shen et al. (2011), with the main difference being that we constrain the shape of the H $\beta$ line to be the same as that of H $\alpha$. Due to the large amount of dust extinction in these reddened quasars, the S/N is considerably higher in H $\alpha$ than in H $\beta$ and we find that we obtain better fits to the region around H $\beta$ by constraining the shape of the emission line in this way.

3.1.1 H $\alpha$ modelling procedure

A power law continuum is fit to the data in the wavelength regions 6000-6350 and 6800-7000 Å, except for VIK J2238-2836, where the power law is instead fit to the region 6800-7800 Å. This continuum is then subtracted from the spectrum and seven different models are fit to the H $\alpha$ line in the window 6350-6800 Å.

We fit the continuum-subtracted H $\alpha$ line with a model consisting of one, two, or three broad Gaussian components. In the case of two Gaussians, this fit is done twice: once constraining both Gaussians to have the same centroid and once where the two components are unconstrained. For the three models with one or two broad components, a separate fit is also tried with the addition of a narrow line region (NLR) template. The NLR template consists of five Gaussians, each constrained to have the same width and velocity offset: one to fit any narrow component of H $\alpha$, and two each to fit the [N ii] $\lambda\lambda$6548,6584 and [S ii] $\lambda\lambda$6717,6731 doublets. The amplitudes of the [S ii] doublet are constrained to be equal and those of [N ii] are fixed at the expected ratio of $1:2.96$. All components of all models are constrained to be non-negative.

Each fit is visually inspected, and models which contain physically unreasonable line widths (\leq 200\hbox{\rm\thinspace km\rm\thinspace s^{-1}}) are rejected. A model with more free parameters is accepted in favour of a simpler model only if there is a greater than 10 per cent reduction in the $\chi^{2}$ per degree of freedom in the fitting window. These ‘best’ fits are shown in Figs. 1 and 12, and have a median reduced $\chi^{2}$ value of 0.978. One object (VIK J2243-3504) is found to require the inclusion of the NLR template in order to best describe the shape of the H $\alpha$ line profile; all the other objects are adequately described by up to three broad Gaussians.

3.1.2 H $\beta$ and [O iii] modelling procedure

The region outside the observed frame wavelengths 1.45-1.84 $\mu\textrm{m}$ is first masked to exclude the portions of the spectrum which are severely affected by atmospheric absorption. The rest frame wavelength region 4435-5535 Å is then fit simultaneously with (i) a power law, (ii) an iron template taken from Boroson & Green (1992), (iii) the profile of the best-fitting H $\alpha$ model, with the normalisation free to vary to fit H $\beta$, and (iv) one of three models for the [O iii] $\lambda\lambda$4960,5008 doublet. The iron template is first convolved with a Gaussian kernel, the width of which is allowed to vary up to the width of the profile of the best-fitting H $\alpha$ model. In all objects, the resulting fits have an iron width consistent with that of the H $\alpha$ emission. We show in Appendix B that our results are robust when comparing the use of different iron templates.

The [O iii] $\lambda$4960 line is constrained to have the same shape as the $\lambda$5008 line, with the total flux in the lines fixed at the theoretical ratio of $1:2.98$ (Storey & Zeippen, 2000). We try fitting zero, one, and two Gaussians to the $\lambda$5008 line, with an extra Gaussian accepted only if it leads to a greater than 10 per cent reduction in the $\chi^{2}$ per degree of freedom in the wavelength range 4700-5100 Å. A fourth [O iii] model with a third Gaussian component is found not to lead to a significant reduction in reduced $\chi^{2}$ in any of the objects in this sample.

The $\lambda$4960 and $\lambda$5008 lines are fit simultaneously, although we only use the $\lambda$5008 line properties when presenting our results. Eleven objects are adequately fit with a single Gaussian describing the $\lambda$5008 line. Eight objects require a second Gaussian component to fully reproduce the $\lambda$5008 line profile, providing evidence for an asymmetric line profile in $\approx 42$ per cent of the population, similar to the fraction found in unobscured quasars by Coatman et al. (2019). These fits are shown in Figs. 1 and 14, and have a median reduced $\chi^{2}$ value of 1.113 across the wavelength range 4700-5100 Å.

Three objects are found to be consistent with no [O iii] emission. We do not make any estimate of the velocity width of the [O iii] lines in these cases. In ULAS J1216-0313 and ULAS J1234+0907, we estimate an upper limit on the equivalent width of the [O iii] emission line by fitting a template, generated by running our fitting routine over a median composite spectrum constructed from the 19 objects in our sample with reliable [O iii] measurements. In one object, ULAS J0144-0014, we find no [O iii] emission, even when constraining the shape of the line with the template. These three objects are shown in Fig. 16.

3.1.3 Black hole masses

For intrinsically luminous, heavily reddened quasars at $z>2$, it is known that host galaxy contamination does not significantly affect the $JHK$ colours (Wethers et al., 2018). Therefore, for the 22 quasars in our high redshift sample, we derive extinctions using the near infra-red photometric data from VISTA and UKIDSS. Where $J$-band photometry is available, we find the $E(B-V)$ required to best fit the observed $J-K$ colour when reddening a quasar template (Maddox et al., 2008; Wethers et al., 2018). If $J$-band photometry is not available, we use the $H-K$ colour instead. Rest frame 5100 Å luminosities are corrected for dust extinction using this $E(B-V)$. We derive black hole masses using the prescription from Vestergaard & Peterson (2006), and use the correction from Coatman et al. (2017) to estimate the full width at half maximum (FWHM) of H $\beta$ from our fit to H $\alpha$:

[TABLE]

[TABLE]

Eddington ratios are estimated assuming $L_{\rm bol}=8\times L_{5100}$, and are given in Table 3, alongside the derived extinctions, dust-corrected luminosities and black hole masses.

We find extinctions in the range $0.5<E(B-V)<2.0$, as expected from our photometric selection criteria. For our faint sample, we find black hole masses in the range 10^{8.4-9.4}\hbox{\rm\thinspace M_{\odot}}, generally smaller than those presented in Banerji et al. (2012); Banerji et al. (2015). These quasars are intrinsically fainter than previously known heavily reddened quasars and do not simply appear fainter due to increased extinction. For our bright sample, we find more massive black holes in the range 10^{9.2-10.5}\hbox{\rm\thinspace M_{\odot}}, placing them among the most massive black holes known at these redshifts (e.g. Bischetti et al., 2017).

3.2 Low redshift sample

Twelve of the 25 spectroscopically confirmed quasars presented in Section 2 are at $z<2$. At these lower redshifts, the observed colours may be contaminated by emission from the host galaxy, especially in bluer bands (e.g. $J$) where (by selection) the quasar is faint. To mitigate against this source of uncertainty, we estimate the extinction by comparing the $H-K$ colour from our SINFONI spectra (which are extracted from smaller apertures and will thus be less affected by host galaxy emission) with the colours predicted from reddening a quasar template (Maddox et al., 2008; Wethers et al., 2018). We find extinctions in the range $1.3\leq E(B-V)\leq 3.0$, and correct the line luminosities for these extinctions before deriving black hole masses. Results are shown in Table 2.

For objects with redshifts in the range $1.3<z<2.0$, we fit a single Gaussian to H $\alpha$ and estimate the black hole mass using the scaling relation from Greene & Ho (2005):

[TABLE]

We find FWHMs in the range 2000-6500$\rm\thinspace km\rm\thinspace s^{-1}$and black hole masses in the range 10^{8.6-9.6}\hbox{\rm\thinspace M_{\odot}}.

For objects with $z<1.3$, we use the scaling relation between the luminosity and width of the Pa $\beta$ emission line and the black hole mass given by method 2 of Kim et al. (2010):

[TABLE]

For objects with $z<0.9$, we fit a single Gaussian to the Pa $\beta$ line. For objects with $0.9<z<1.3$ where we have coverage of neither H $\alpha$ nor Pa $\beta$, we fit kinematically tied Gaussians simultaneously to Pa $\gamma$, Pa $\delta$, Pa $\epsilon$, He i $\lambda 10830$, and O i $\lambda 11287$, with the amplitude of each line free to vary. To convert $L_{\rm{Pa}\,\gamma}$ to $L_{\rm{Pa}\,\beta}$, we assume a flux ratio of Pa $\gamma:\rm{Pa}\,\beta=1:1.80$ (Storey & Hummer, 1995; Dopita & Sutherland, 2003). We find Paschen line FWHMs in the range 1000-4500$\rm\thinspace km\rm\thinspace s^{-1}$and black hole masses in the range 10^{7.5-9.1}\hbox{\rm\thinspace M_{\odot}}.

These lower redshift quasars are therefore typically more obscured and have slightly smaller black hole masses than heavily reddened quasars at $z>2$, such as the objects from Banerji et al. (2012) and the high redshift sample in this paper.

4 Emission line properties

In this section we explore the [O iii] $\lambda$5008 emission line properties in our high redshift sample (Section 2.3). We use the best-fitting [O iii] model as described in Section 3.1.2 to measure the shape and equivalent width of the line. For those objects which require multiple Gaussian components to best model the emission line, we do not ascribe any significance to individual components, but instead derive properties for the total line profile.

We compare the objects in our high redshift sample to a large sample of unobscured quasars in the same redshift and luminosity range. This comparison sample contains 111 quasars with $2.0<z<2.6$ and near infra-red spectra taken from the catalogue of Coatman et al. (2019, see Fig. 2), which are flagged in the catalogue as having robust measurements of the [O iii] emission line properties. The selection of quasars in this catalogue is described in detail by Coatman et al. (2017). In the context of this work, the comparison sample essentially provides an unbiased representation of the optical SEDs of unobscured quasars, matched in redshift and luminosity to our high redshift reddened quasar sample. It is worth noting that using a larger comparison sample taken from the full catalogue redshift range of $1.5<z<4.0$ does not alter any of the results in our subsequent analysis.

To help visualise the quality of data in the two samples, in Fig. 3 we show stacks of our spectra in the region around the [O iii] emission line, normalised at 5100 Å. For the purposes of this stack only, we take a subsample of 22 objects from our comparison sample which are the nearest neighbours to each of our 22 reddened quasars in the redshift-luminosity space shown in Fig. 2. We apply no correction for dust obscuration before stacking, and the slope of the continuum across the narrow wavelength range shown is clearly redder in the reddened sample. Comparing the two stacks, it can be seen that (a) the reddened objects generally have lower S/N in this part of the spectrum, (b) the two samples have very similar [O iii] line widths, and (c) the [O iii] emission in the reddened sample is perhaps slightly stronger relative to the continuum. We will explore the shapes and strengths of the [O iii] lines in our two samples further in Sections 4.1, 4.2 and 4.3.

Uncertainties on derived parameters are estimated using a Monte Carlo approach. One hundred realisations of each spectrum are generated by smoothing the observed quasar spectrum with a 5-pixel inverse-variance weighted top-hat filter, and then adding ‘noise’ drawn from a Gaussian distribution with dispersion equal to the spectrum flux uncertainty at each wavelength. Quoted parameter uncertainties are 1.48 times the median absolute deviation of the parameter distribution in the ensemble of simulations.

We compare the samples using the two-sample Kolmogorov–Smirnov (KS) test (Peacock, 1983), which tests the null hypothesis that two samples are drawn from the same underlying probability distribution (i.e., the same population). We also use the two-sample Anderson–Darling (AD) test (Anderson & Darling, 1952; Darling, 1957), which tests the same null hypothesis, but is more sensitive to differences in the tails of the distributions.

We begin by comparing the strengths of the broad H $\alpha$ emission in our samples, as measured by the equivalent width. Our heavily reddened quasars were selected by virtue of their red colours, leading to a selection bias towards objects which are bright in the $K$-band, i.e. those objects which might have stronger H $\alpha$ emission in the redshift range $2.0<z<2.6$. As shown in Fig. 4, we find no systematic difference in the equivalent width of the broad line emission between the two samples, suggesting that the amount of ionised gas is not significantly different in the reddened and unobscured samples.

We note that, if one adopts a theoretical value for the unreddened H $\alpha$ to H $\beta$ ratio, the observed ratio can be used to estimate the extinction in front of the broad line emitting region. However, in addition to the statistical uncertainty of $\sim 0.1$ magnitudes arising from the low S/N in the H $\beta$ lines of our heavily reddened quasars, further uncertainty arises due to the unknown temperatures and densities of the gas producing the Balmer emission, i.e. the assumption of a fixed intrinsic flux ratio (Korista & Goad, 2004). Adopting an intrinsic ratio of $\rm{H}\,\alpha:\rm{H}\,\beta=2.86:1$ (Dopita & Sutherland, 2003), we find that the extinctions estimated in this way are consistent with the $E(B-V)$s we derive from photometry in Section 3.1. However, at a fixed value of the $E(B-V)$ derived from photometry, the scatter in the extinction we measure from the Balmer lines is $\sim 0.5$ magnitudes. We therefore use $E(B-V)$s from photometry in the rest of this paper.

4.1 The shape of the [O iii] lines

As described in Section 3.1.2, we have information about the shape of the [O iii] lines in 19 of our 22 high redshift reddened quasars. For these 19 objects, we first calculate the 80 per cent velocity width ($w_{80}$) of the $\lambda$5008 line profile. We take the best-fitting [O iii] model and integrate to find the total flux under the line profile. The velocity width is then defined to be $w_{80}=v_{90}-v_{10}$, where $v_{90}$ and $v_{10}$ are the respective velocities at which 90 and 10 per cent of the cumulative flux in the line profile is found. We correct for instrumental broadening of \sigma=200\hbox{\rm\thinspace km\rm\thinspace s^{-1}}, noting that all of our [O iii] lines are well resolved with the smallest velocity width in our sample being w_{80,\rm{min}}=600\hbox{\rm\thinspace km\rm\thinspace s^{-1}} prior to applying the correction.

The [O iii] velocity widths in our reddened sample and in the comparison sample of unobscured quasars are shown in Fig. 5. We find a weak correlation between $w_{80}$ and continuum luminosity, as is also found in unobscured quasars at the same redshifts and luminosities (Netzer et al., 2004; Shen, 2016; Coatman et al., 2019). There is no such correlation between $w_{80}$ and either the black hole mass, the Eddington ratio or the $E(B-V)$. Using the KS and AD tests, we find no evidence to reject the null hypothesis that the two samples are drawn from the same underlying distribution of velocity widths. As the velocity widths in both samples are significantly larger than the velocity dispersion expected from virialised gas motions, it is reasonable to assume that most, it not all, of our velocity width measurements are dominated by emission from outflowing gas. Our results therefore suggest that the ionised gas outflows in the host galaxies of our reddened quasars are, on average, no faster or slower than those in the unobscured population.

In Section 3.1.2, we found that around 42 per cent of the quasars in both the reddened and unobscured samples show evidence for asymmetric [O iii] emission, in that they are not adequately described by a single Gaussian line profile, and require a second Gaussian component in the fit to [O iii]. For these objects, we quantify the asymmetry of the $\lambda$5008 line using $(v_{90}-v_{50})/(v_{50}-v_{10})$, and find no significant difference between the two samples in terms of the distribution of the amount of $\lambda$5008 line asymmetry.

4.2 The equivalent width of the [O iii] lines

Due to the uncertainty in the location and physical scale of the [O iii] emitting material and also the uncertainty in the location of the obscuring dust in heavily reddened quasars, the amount by which the [O iii] emission in our sample is attenuated is unknown. We can assume that the [O iii] attenuation is at least zero, but no more than the attenuation of the continuum: the obscuring dust might be located entirely in front of the [O iii] emitting region, or behind the [O iii] emitting region but in front of the continuum source. This unknown factor in the geometry of these objects leads to an uncertainty when deriving any measure of the intrinsic strength of the [O iii] emission.

With this uncertainty in mind, we derive the equivalent width (EW) of the [O iii] $\lambda$5008 line under two different assumptions: that the emission line is reddened by the same amount as the continuum, and that the continuum is reddened but that the emission line itself is completely unobscured. The first case will produce the largest possible EW, and the second will give the smallest possible EW, as de-reddening the continuum while leaving the emission line unchanged increases the flux in the continuum and hence lowers the EW. Any amount of ‘intermediate’ de-reddening of the [O iii] emission will produce an EW somewhere between these two extremes.

In both cases, the cumulative distribution function (CDF) of the [O iii] EWs in our sample is formed and shown in Fig. 6. For the case where the [O iii] emission is unobscured and so does not need to be corrected for extinction, we find that we can reject the null hypothesis with $p<0.001$, meaning that such a distribution of [O iii] EWs is inconsistent with having been drawn from the same population as the unobscured quasars in Coatman et al. (2019).

However, we are unable to reject the null hypothesis ($p>0.3$) in the case that the obscuring dust is located completely in front of the [O iii] emitting gas, meaning that our results are consistent with a scenario in which the intrinsic distribution of [O iii] EWs in our reddened sample is the same as that of the unobscured population, and where the amount of obscuration in our objects along the line of sight towards the [O iii] emitting gas is is the same as the amount of obscuration along the line of sight to the continuum source.

4.3 Degeneracy between the amount of ionised gas and the location of the reddening dust

When deriving any measure of the strength of the [O iii] emission line, such as the equivalent width, there is a degeneracy between the intrinsic amount of line emission and the amount of obscuration along the line of sight towards the emitting material. The measured equivalent width of the emission line will always depend on both of these factors. For the heavily dust-reddened objects presented in this paper, we have no independent constraints on the location of the obscuring material - more specifically, the dust could lie inside or outside the [O iii] emitting region, or could be co-located with it. Due to this degeneracy, we must make an assumption in order to draw any further conclusion about the relative strength of the ionised gas emission in reddened quasars compared to the unobscured population in the same redshift and luminosity range.

As the velocity widths of the [O iii] lines, the continuum luminosities, and the H $\alpha$ equivalent widths in our sample are consistent with having been drawn from the same distributions as those of the unobscured population, it is reasonable to assume that the amount of ionised gas in our reddened quasars is the same as in the unobscured population and so the intrinsic [O iii] emission is not significantly stronger or weaker in the reddened population.

If we make this assumption, i.e. that there is no difference in the distribution of the amount of ionised gas in the two populations and hence that the two populations have the same intrinsic [O iii] EW distribution, then we can fit the observed CDF of our reddened sample to that of the unobscured population, as shown in Fig. 7. We find that the amount of extinction required on average in the [O iii] lines of our reddened objects to best fit the EW distribution is 0.92 times the $E(B-V)$ derived from the near infra-red photometry in Section 3.1.3: on average, just 8 per cent of the continuum extinction arises inside the [O iii] emitting region. However, it is worth reemphasising the result in Section 4.2, in that, under the assumption that the two populations have the same [O iii]  EW distribution, our reddened quasars are not inconsistent with having the obscuring material located entirely outside the [O iii] emitting region.

We note that the two most luminous quasars in our reddened sample (ULAS J1216-0313 and ULAS J1234+0907) are the two most highly reddened objects and also have very weak [O iii] emission. Even if we correct the observed [O iii] for the largest possible extinction, the EWs in these two objects are still less than half the median EW in the rest of the reddened sample, suggesting that the weakness of the observed [O iii] emission in these objects cannot be explained purely by obscuration, and that the [O iii] emission is intrinsically weak in these objects. This result is very similar to the fraction ($\sim 10$ per cent) of unobscured quasars at these luminosities and redshifts which have very weak [O iii] emission (i.e. EW$<1$Å; Coatman et al., 2019), and we suggest that it is consistent with a scenario in which the most luminous quasars show weaker emission from the narrow line region (Netzer et al., 2004), either by overionising or by physically removing the emitting gas.

5 Discussion

5.1 Ionised gas kinematics

In Section 4.1, we find no significant difference between the kinematics of the [O iii] emitting gas in the heavily reddened and unobscured quasar populations at $z>2$, as traced by either the velocity width or the asymmetry of the [O iii] line profile, although we do note that there would need to be a large dissimilarity between the two populations in order to find a statistically significant difference between our samples, due to the small sample sizes. Thus, while the kinematics of the [O iii] emission are not observed to be significantly different, we cannot rule out the possibility of smaller differences in the ionised gas emission properties of unobscured and heavily reddened quasars. We also show that the [O iii] velocity width tends to increase with increasing AGN luminosity in both populations, a trend that has also been seen in lower redshift AGN (Zakamska & Greene, 2014; Harrison et al., 2016).

Wu et al. (2018) report the detection of [O iii] $\lambda\lambda$4960,5008 in W1136+4236, a WISE-selected Hot Dust-Obscured Galaxy (HotDOG). W1136+4263 is at $z=2.41$, with a dust-corrected luminosity of L_{5100}\approx 10^{46.5}$$\rm\thinspace erg\rm\thinspace s^{-1}\,and an extinction of $E(B-V)=2.5$: while it is more reddened than the reddened quasars in our sample, it lies in the same luminosity and redshift range. This HotDOG has an [O iii] velocity width of w_{80}\approx 2500\hbox{\rm\thinspace km\rm\thinspace s^{-1}}, placing it within the range spanned by our data in Fig. 5. It is hard to draw conclusions from this single object, which is the least obscured HotDOG studied by Wu et al., but there is currently no evidence to suggest that the ionised gas kinematics in HotDOGs are any more extreme than those in our heavily reddened quasar population. Similarly, the red quasars presented by Urrutia et al. (2012) with L_{\rm bol}\approx 10^{46}$$\rm\thinspace erg\rm\thinspace s^{-1}\,show [O iii] emission with the width of the broad component up to 1600$\rm\thinspace km\rm\thinspace s^{-1}$(Brusa et al., 2015), which coincides with the area of parameter space populated by our heavily reddened quasars.

Comparing to type 2 quasars, we note that the ionised gas kinematics in our reddened quasars are similar to the [O iii] velocity widths in the sample of type 2 quasars at $z<0.8$ studied by Zakamska & Greene (2014), who find $w_{80}=280$-3000$\rm\thinspace km\rm\thinspace s^{-1}$. These quasars have L_{[\text{O\,{iii}}]}=10^{42-43.5}$$\rm\thinspace erg\rm\thinspace s^{-1}\,. However, the broad line and continuum emitting regions in these objects are obscured by dust on nuclear scales which is unlikely to be obscuring the [O iii] emitting region, and so it is hard to directly compare the luminosities of these quasars with the objects in our sample.

Our results complement those of Harrison et al. (2016), who examined the [O iii] properties of an X-ray selected sample of optically unobscured AGN at $1.1<z<1.7$. The AGN in their sample have [O iii] $w_{80}$ in the range 200-1000$\rm\thinspace km\rm\thinspace s^{-1}$, L_{[\text{O\,{iii}}]}\approx 10^{42}\hbox{\rm\thinspace erg\rm\thinspace s^{-1},} and L_{X}\approx 10^{44}$$\rm\thinspace erg\rm\thinspace s^{-1}\,. They are therefore around an order of magnitude fainter than our heavily reddened quasars (Lansbury et al. in prep.) and also appear to have narrower [O iii] line profiles in comparison to our quasars. Harrison et al. (2016) found no significant difference between the ionised gas kinematics in X-ray obscured and X-ray unobscured objects. Taken with our results, this supports a scenario where the most powerful ionised gas outflows in massive galaxies are driven by the luminosity of the central quasar, independently of the column density of any obscuring material along the line of sight to the observer.

By contrast, DiPompeo et al. (2018) do find a difference in the [O iii] properties of red and blue colour-selected AGN, with the redder AGN having broader, more blueshifted [O iii] emission on average. However, the DiPompeo et al. AGN are significantly less luminous and at much lower redshifts than the quasars presented in this work. The two results are therefore not inconsistent and could indicate different redshift evolution in the outflow properties of reddened and unobscured AGN.

We note that the heavily reddened quasars in our sample do not show the same extreme [O iii] kinematics as the four rare (i.e. out of 97 ERQs; Hamann et al. 2017) objects presented in Zakamska et al. (2016), which are in the same luminosity and redshift range as our heavily reddened quasars. However, the Zakamska et al. quasars were also selected to have high equivalent width C iv emission, and it is conceivable that the extreme [O iii] kinematics seen in that sample could be linked to this selection. Alternatively, Coatman et al. (2019) find similarly rare (18/354) unobscured quasars with ionised gas kinematics that are almost as extreme, and so it could be that $\lesssim 5$ per cent of the most luminous quasars at $z>2$ have ionised gas outflows of up to 5000$$\rm\thinspace km\rm\thinspace s^{-1}, irrespective of reddening, and that if we were to increase our sample of luminous heavily dust-reddened quasars by a factor of $\sim 5$ we might also observe such extreme outflows.

If the objects in our reddened sample do indeed represent the ‘blow-out’ phase for luminous quasars at redshifts $z\approx 2$, then, as the signatures of fast outflows are also seen in the unobscured population at the same redshift, our results suggest that the outflowing gas persists after the obscuring dust has been cleared from the line of sight. In particular, the asymmetries and large velocity widths seen in the [O iii] emission in our comparison sample of unobscured quasars would be due to outflows which have already cleared any obscuring material away from the line of sight to the continuum source.

5.2 Obscuring dust

In Sections 4.2 and 4.3, we found that the obscuring dust in our high redshift sample of heavily reddened quasars was most likely to be located almost entirely outside the [O iii] emitting region, suggesting that the obscuration in heavily reddened quasars arises on kiloparsec scales. Using data from ALMA, Banerji et al. (2018) find large (i.e. galaxy-scale) dust emitting regions in at least some heavily reddened quasars, consistent with this result. Deeper observations providing spatially resolved information about either the location of the dust or the location of the ionised gas are necessary to place firmer constraints on the physical scales and covering factor of the dust relative to the ionised gas in our sample.

We note that our results could be interpreted with a model whereby our reddened quasars are in fact unobscured quasars which are viewed along a line of sight that is skimming the plane of the host galaxy. However, we disfavour such a model as it would (a) require a significant disc component to the host galaxy, which would not be consistent with a merger-driven AGN triggering event, and (b) require the [O iii] emitting gas to be located very close in to the centre of the galaxy, in contrast with local AGN (see Section 1).

Our results are consistent with a model of galaxy formation and evolution where outflows persist long after the obscuring dust has cleared from the line of sight. There are also examples of unobscured, UV-luminous quasars at high redshifts with significant amounts of cool dust inferred from their FIR and sub-millimeter continuum emission (Harris et al., 2016; Pitchford et al., 2016). The general consensus is that this cool dust is being heated by star formation in the quasar host galaxy rather than the quasar itself (although see Symeonidis 2017 for an alternative explanation). Thus, star formation in luminous quasar hosts can also persist after the dust has been cleared from the line of sight. Using a population of SDSS quasars, Maddox et al. (2017) find that strong FIR emission (and therefore presumably star formation) is associated with strong nuclear outflows as traced by large C iv blueshifts. Taken together with the results of Coatman et al. (2019), where a correlation between the blueshift of [O iii] and the blueshift of the C iv emission line is seen, this implies that FIR-luminous quasars also exhibit stronger galaxy-scale outflows compared to their FIR-faint counterparts. The FIR bright quasars studied by Maddox et al. (2017) are not reddened, but outflows and starburst activity are nevertheless being seen in this population. There is likely to be a considerable overlap between star formation, outflows, and the presence and clearing of dust in massive (active) galaxies, and obtaining observations at multiple wavelengths of the same objects in order to quantify these processes is the only way to move to a complete understanding of the interplay between SMBH accretion, outflows, and star formation, and their impact on galaxies at $z>2$.

6 Summary

From 26 candidates, selected by virtue of their red colours in the near infra-red, we have spectroscopically confirmed 25 new heavily dust-reddened broad line quasars at $0.7<z<2.6$, which extend the bright and faint ends of the luminosity function presented in Banerji et al. (2015).

By combining with the objects presented in Banerji et al. (2012), we construct a sample of 22 heavily reddened quasars with extinctions $0.5<E(B-V)<2.0$, luminosities $45<{\rm Log}(L_{5100}/\rm\thinspace erg\rm\thinspace s^{-1})<48$, redshifts $2.0<z<2.6$, and spectra covering the H $\alpha$, H $\beta$, and [O iii] emission lines. Using the [O iii] line, we derive parameters describing the ionised gas kinematics and compare to a large sample of 111 unobscured quasars in the same redshift and luminosity range. This is the first time that the ionised gas kinematics in the reddened and unobscured quasar populations have been systematically compared at redshifts $z>2$. Our main results are:

(i)  There is no significant difference between the two populations in the velocity widths or asymmetries of their [O iii] lines, suggesting that ionised gas outflows are moving, on average, no faster or slower in the reddened population than in the unobscured population.

(ii)  As our samples are drawn from the same luminosity and redshift range, this suggests that the outflow driving mechanisms are not significantly different in the reddened and unobscured populations.

(iii)  The [O iii] equivalent width distributions in the two populations are consistent with a scenario in which the reddened quasars have the same strength of ionised gas emission as the unobscured quasars, and the [O iii] emission in the reddened objects is subject to the same amount of extinction as the continuum. In other words, the obscuring dust is most probably located outside the [O iii] emitting region on kpc-scale distances from the central SMBH.

We suggest that our results are best explained by a model in which, following a starburst episode, quasar-driven ionised gas outflows persist for some time after any obscuring dust has been cleared from the line of sight. Such an episode of star formation would most likely be triggered by a major galaxy merger, which would also provide fuel for black hole accretion and hence the trigger for quasar activity.

Acknowledgements

We thank George Lansbury and Richard McMahon for useful discussions. An anonymous referee provided a comprehensive review of the paper for which we are grateful. MJT thanks the Science and Technology Facilities Council (STFC) for the award of a studentship. MB acknowledges funding from the Royal Society via a University Research Fellowship. PCH acknowledges funding from STFC via the Institute of Astronomy, Cambridge, Consolidated Grant. CP thanks the Alexander von Humboldt Foundation for the granting of a Bessel Research Award held at MPA. CP is also grateful to the ESO and the DFG cluster of excellence ‘Origin and Structure of the Universe’ for support.

This work is based on observations collected at the European Southern Observatory under ESO programme 095.A-0094(A), and made use of Astropy (Astropy Collaboration et al., 2013; Price-Whelan et al., 2018) and Matplotlib (Hunter, 2007).

Appendix A Spectra and Emission Line Fits

We show the new SINFONI spectra for our 26 targets in Figure 8, and our fits to H $\alpha$, H $\beta$ and [O iii] in our ‘high redshift sample’ in Figures 1, 12, 14 and 16.

Appendix B Iron Subtraction

We test the sensitivity of our results to the assumed shape of the iron template by replacing the iron template of Boroson & Green (1992, BG92) with that of Kovačević et al. (2010, KPD10) and re-deriving the [O iii] line properties with an otherwise identical line fitting routine. The iron prescription of KPD10 is much more flexible, in that it allows the ratios of different groups of iron lines to vary, corresponding to an increase in the number of free parameters in the fit and enabling a more accurate description of the iron in objects where iron line ratios are different from I Zw 1, the NLS1 galaxy used to derive the template of BG92.

We compare the main results of our work under the two different iron prescriptions, as shown in Fig. 17. We visually inspect individual objects, such as VHS J1117-1528, which has w_{80}=2330,3130$$\rm\thinspace km\rm\thinspace s^{-1} and $\rm{EW}=22.7,33.8$ Å when fit with BG92 and KPD10 respectively, making this the object from our sample with the biggest difference in the derived [O iii] velocity widths. The fit in this quasar (Fig. 18) is sensitive to the form of the iron prescription due to the fact that the [O iii] emission line is both broad (w_{80}>2000$$\rm\thinspace km\rm\thinspace s^{-1}) and weak ($\rm{EW}<35$ Å), and we caution other authors to be careful when deriving line properties in similarly extreme objects.

However, the population statistics in our sample do not change significantly when using a different iron prescription to fit our data; in particular, the results described in Section 4 are the same when derived with either prescription. Therefore, while we are cautious about the iron contribution in individual objects, we are confident that the precise form of the iron template used does not affect the conclusions of this paper.