New Constraints on Quasar Evolution: Broad Line Velocity Shifts over $1.5\lesssim z\lesssim 7.5 $

TL;DR

This study analyzes UV spectra of 394 quasars across redshifts 1.5 to 7.5, revealing that C IV line blueshifts remain stable at lower redshifts but increase significantly at z > 6, suggesting evolving quasar wind properties.

Contribution

Introduces a model-independent method for measuring quasar emission line shifts and provides new insights into quasar wind evolution over cosmic time.

Findings

C IV blueshift does not significantly evolve from z=1.5 to 6

C IV blueshift increases by a factor of ~2.5 at z > 6

Method yields unbiased velocity shifts insensitive to spectral resolution

Abstract

We present the results of a model-independent investigation of the rest-frame UV spectra from a comprehensive sample of $394$ quasars in the redshift range $1.5\leq z \leq 7.5$. We fit the main Broad Emission Lines (BELs) in the rest-frame range $1280 \text{\AA} \leq \lambda \leq 3000 \text{\AA}$ (O I, C II, Si IV, C III, C IV and Mg II) with a lightly-supervised spline fitting technique. Redshifts are derived from the peaks of each fitted BEL and used to compute relative velocity shifts. We show that our method gives unbiased velocity shifts and is insensitive to spectral resolution and instrumental parameters. It is found that the average blueshift of the \cfour\, line with respect to several low-ionisation lines in luminosity-matched samples does not significantly evolve over $1.5\gtrsim z\gtrsim6$. However, the average blueshift increases significantly by a factor $\sim 2.5$ at $z\gtrsim 6$. We propose that this redshift evolution can be explained by \cfour\, winds launched perpendicularly to an accretion disk with increased torus opacity at high-redshift, coupled with a potential orientation-driven selection bias. Our results open new exciting avenues of investigation into young quasars in the reionisation epoch.

Figures

Tables

Table 1. Table 1: The different quasar samples with their redshift and luminosity ranges, resolution, median SNR and number of objects. (a) We use here the subset of DR12Q that have both tabulated C  IV   and Mg  II -derived redshift. (b) Throughout the text, SDSS refers to a subsample of the DR7Q and DR12Q overlapping quasars, luminosity-matched to the z6 sample, described in Section 2.1

Name	$\Delta z$	$\Delta L_{\lambda ,1350}$	R	$\overline{\text{SNR}}$	$N$
DR12Qa	$1.5-2.4$	40.6-46.0	1500-2500	3.4	99232
SDSSb	$1.5-4.5$	43.3-44.4	1500-2500	14	108
XQ100	$3.5-4.5$	44.0-45.2	6000-9000	30	100
GGG	$4.4-5.5$	42.7-44.3	800-900	20	163
z6	$5.4-6.3$	43.4-44.3	9000-11000	90	11
z7	$6.4-7.5$	42.8-43.8	$\sim 6000$	6	12

Table 2. Table 2: Median ( Δ ​ v ¯ ¯ Δ 𝑣 \overline{\Delta v} ), mean ( ⟨ Δ ​ v ⟩ delimited-⟨⟩ Δ 𝑣 \langle\Delta v\rangle ) and standard deviation ( σ Δ ​ v subscript 𝜎 Δ 𝑣 \sigma_{\Delta v} ) of the velocity error Δ ​ v Δ 𝑣 \Delta v on of the recovered peak of each BEL with respect to the original value derived from XQ100 spectra. The columns one to four (‘Rebinned’) give the values when the resolution is lowered to SDSS-like levels by rebinning the data whereas the next three (‘SDSS’) show the comparison between redshifts retrieved from XShooter and SDSS spectra of the same quasars. All values are provided in km s -1 .

	Rebinned			SDSS
Line	$\overline{\Delta v}$	$⟨\Delta v⟩$	${\sigma }_{\Delta v}$	$\overline{\Delta v}$	$⟨\Delta v⟩$	${\sigma }_{\Delta v}$
O I	$-56$	$-44$	$250$	$-73$	$-38$	$349$
C II	$90$	$-6$	$538$	$-25$	$-6$	$179$
Si IV	$-96$	$-33$	$451$	$-53$	$-11$	$407$
C IV	$-30$	$-27$	$537$	$-30$	$-16$	$131$
C III]	$-57$	$-30$	$474$	$36$	$19$	$236$
Mg II	$25$	$-0$	$166$	$-35$	$-19$	$209$

Table 3. Table 3: BEL-based redshifts for all our objects. SDSS quasars are designated by their plate-mjd-fiber identification to make their search easier in the SDSS database. The full table is published as an online supplementary material.

Quasar	$z_{\text{O}\text{I}}$	$z_{\text{C}\text{II}}$	$z_{\text{Si}\text{IV}}$	$z_{\text{C}\text{IV}}$	$z_{\text{C}\text{III}]}$	$z_{\text{Mg}\text{II}}$
SDSS Quasars
5379-55986-0904	3.87	3.87	3.87	3.81	3.84
4481-55630-0178	2.52		2.51	2.50	2.49	2.52
4193-55476-0376	2.28		2.26	2.26	2.25	2.27
$⋮$		$⋮$		$⋮$		$⋮$
XQ100 Quasars
HB89 0000-263	4.13		4.11	4.10	4.12
BR J0006-6208	4.46			4.45	4.44
BR J0030-5159	4.18		4.17	4.18	4.17
$⋮$		$⋮$		$⋮$		$⋮$
GGG Quasars
SDSSJ0011+1446	4.96	4.97	4.95	4.94
SDSSJ0822+1604			4.47	4.50
SDSSJ1043+6506			4.48	4.46
$⋮$		$⋮$		$⋮$		$⋮$
z6 Quasars
J0148+0600	5.97	5.96	5.93	5.91	5.93	5.98
J0836+0054	5.80	5.80	5.81	5.78	5.76	5.77
J0927+2001	5.74	5.76	5.73	5.72	5.77	5.76
J1030+0524	6.30	6.27	6.30	6.28	6.31	6.30
J1306+0356	6.03	6.02	6.02	6.00	6.01	6.02
J1319+0950	6.10	6.11	6.17	6.03	6.14	6.12
J0100+2802	6.31		6.25	6.24	6.30
J0818+1722					5.95
J1509-1749	6.08	6.11	6.11	6.10	6.10	6.12
J1044-0125	5.78	5.80	5.77	5.76	5.79	5.78
J0231-0728	5.42	5.42	5.43	5.42	5.42	5.42
z7 Quasars
J2348-3054	6.91	6.89	6.87	6.85	6.89	6.89
P231-20	6.57		6.68	6.47	6.54	6.60
P167-13	6.50	6.49	6.41	6.41	6.48	6.51
P036+03			6.45	6.42	6.50	6.53
J0305-3150	6.56		6.57	6.57		6.61
P183+05		6.48	6.40	6.35		6.42
J1120+0641		7.09	7.04	7.01	7.06	7.09
J1342+0928	7.56	7.55	7.41	7.35		7.52
P247+24	6.45	6.46	6.44	6.42		6.48
P338+29		6.64	6.70	6.63	6.64
P323+12	6.61	6.58	6.59	6.58	6.59	6.58
J0109-3047	6.77	6.75	6.76	6.69		6.74

Equations

$$F(\lambda )=F_0{\left(\frac{\lambda }{2500\mathring{A}}\right)}^{\alpha }$$

Full text

††pubyear: 2018††pagerange: New Constraints on Quasar Evolution: Broad Line Velocity Shifts over $1.5\lesssim z\lesssim 7.5$–C

1 Introduction

The number of known quasars in the early Universe has recently increased sharply, mainly due to the availability of deep wide-field surveys in the optical and the near-infrared. At the time of writing, we count more than $170$ quasars at $z>6$, and $5$ at $z>7$. This progress is due to the concerted efforts of different teams selecting and confirming candidates from diverse surveys such as the Sloan Digital Sky Survey (SDSS, Jiang et al., 2016), Panoramic Survey Telescope and Rapid Response System data (Pan-STARSS, Kaiser et al., 2010; Bañados et al., 2016), the Dark Energy Survey - Visible and Infrared Survey Telescope for Astronomy (VISTA) Hemisphere Survey (DES-VHS, Reed et al., 2015), the VISTA Kilo-Degree Infrared Galaxy Survey (VIKING, Venemans et al., 2013), the VLT ATLAS survey (Carnall et al., 2015), the UKIRT Infrared Deep Sky Survey (UKIDSS, Venemans et al., 2007; Mortlock et al., 2009, 2011), the DESI Legacy imaging Surveys (DELS) (Dey et al., 2019), the Wide-field Infrared Survey (WISE) survey, and, increasingly, from overlaps of some thereof (Bañados et al., 2018; Wang et al., 2018a, b; Yang et al., 2018; Fan et al., 2019). The recent discovery of a gravitationally lensed quasar may suggest that we are missing further high-redshift quasars (Fan et al., 2019; Pacucci & Loeb, 2019), and therefore it is possible their numbers will swell in the coming years.

This surge has enabled a large number of investigations into a variety of topics. Rising numbers of high-redshift quasars offer the possibility to study their central supermassive black holes. The high masses of the high-redshift black holes powering these quasars have proven a particular challenge to explain as they imply either high seed masses and/or continuous accretion at nearly the Eddington rate for most of their young life (e.g. Willott et al., 2010; Bañados et al., 2018). The study of the co-evolution of galaxies and quasars has also benefited from larger samples of younger quasars (e.g. Venemans et al., 2012; Decarli et al., 2017). Determination of the quasar luminosity function at $z\gtrsim 5.5$ has constrained the contribution of quasars to the cosmic reionisation of hydrogen (e.g. Onoue et al., 2017; Kulkarni et al., 2018). Used as the most distant beacons of light in the Universe, quasars are also an invaluable tool to study the intervening material in absorption. Metal absorption studies have charted the cosmic densities of metals up to the first billion years (e.g. Becker et al., 2009; Ryan-Weber et al., 2009; D’Odorico et al., 2013; Bosman et al., 2017; Codoreanu et al., 2018; Meyer et al., 2019). The scatter in the Lyman-$\alpha$ forest at the tail end of reionisation has highlighted the patchy nature of this process (Bosman et al., 2018; Eilers et al., 2018a). Finally, reduced absorption of the intrinsic Lyman-$\alpha$ emission by neutral gas in the so-called quasar ‘near-zone’ has been explored as a promising test of the evolving neutral fraction of hydrogen at $z\gtrsim 6$ (Mortlock et al., 2011; Keating et al., 2015; Eilers et al., 2017; Bañados et al., 2018).

However, absorption studies rely on the accurate reconstruction of the quasar continuum. It is common practice to extrapolate the Lyman-$\alpha$ emission and the blueward continuum (both heavily absorbed in the reionisation era) from the redward continuum. The most common techniques are based on the fitting of spectral Principal Component Analysis (PCA) eigenvectors learnt from SDSS quasars (e.g. Suzuki, 2006; Lee et al., 2012; Pâris et al., 2011; Davies et al., 2018). Recent studies have debated whether or not these SDSS-PCA templates are suitable for early quasars (Mortlock et al., 2011; Bosman & Becker, 2015; Bañados et al., 2018; Venemans et al., 2018; Davies et al., 2018). A crucial question is hence whether the spectra of quasars evolve with redshift. On the one hand, absorption studies are based on the a priori characterisation of the intrinsic continuum at low-redshift and the assumption that it does not evolve. On the other hand, the intrinsic evolution of the spectra of quasars would provide a key observation to constrain the evolution of Active Galactic Nuclei (AGN), their physics and their co-evolution with their host galaxies.

A prime candidate for testing the evolution of quasar spectra is the broad C IV  $\lambda 1549$ Å  emission line. The equivalent width of the C IV  line strongly anti-correlates with the continuum luminosity - the so-called Baldwin effect (Baldwin, 1977). The C IV  line is also known to present a variety of asymmetrical shapes and blueshifts, measured with respect to the quasar’ systemic redshift, are found to be correlated with its equivalent width (e.g. Richards et al., 2011, and references therein). The origin of these peculiarities is often attributed to radiation-driven outflows or poorly understood Broad Line Region physics (BLR). The bulk population of low-redshift quasars (e.g. SDSS, Shen et al., 2008; Pâris et al., 2017) present on average a mildly blueshifted C IV  line ($\sim 800$ km s*-1*). Therefore, the extreme C IV  blueshifts ($\gtrsim 3000$ km s*-1*) reported in the three first $z>7$ quasar (Mortlock et al., 2011; Bañados et al., 2018; Wang et al., 2018b) raised interesting questions about the possible evolution of quasar spectra (and the properties of the C IV  line in particular).

It remains to be seen, though, if the objects detected so far at high-redshift are representative of the whole population of early quasars. Recent observational studies of the C IV  line blueshift have provided contrasting evidence for such evolution (Mazzucchelli et al., 2017; Shen et al., 2019). Different methods for measuring the C IV  blueshift have complicated comparison across authors and redshifts (e.g. Coatman et al., 2016). The determination of the systemic redshift against which C IV  velocity shifts are measured has variously been done via the cross-correlation of a quasar template (Hewett & Wild, 2010; Richards et al., 2011), Independent Component Analysis (ICA) decomposition (Allen & Hewett, in prep., used by e.g. Coatman et al. 2016; Reed et al. 2019), fitting PCA eigenvectors (e.g. Mortlock et al., 2011), or using a reliable quasar Broad Emission Line (BEL) such as Mg II  (e.g. Mazzucchelli et al., 2017; Shen et al., 2019; Reed et al., 2019). Similarly, the location of the C IV  line has been variously determined via the line centroid (Coatman et al., 2016; Reed et al., 2019), the flux-weighted central wavelength (e.g. Mortlock et al., 2011) the peak of a multi-Gaussian (or Gauss-Hermite) profile (e.g. De Rosa et al., 2014; Shen et al., 2016, 2019) or the maximum of a PCA template (e.g. Pâris et al., 2011, 2017). Furthermore, quasar BELs can be broken down into broad and narrow components, which appear to shift differently (Greig et al., 2017), and absorption in the red or blue wing of the emission biases the measurement of the intrinsic shift. Some of these techniques may even hamper the opportunity to probe the possible evolution of quasar spectra as they explicitly rely on templates extracted solely from low-redshift objects. Such factors cast a shadow of uncertainty on these velocity shifts.

In this paper, we aim to address the question afresh and study the evolution of quasar spectra across cosmic time. In order to achieve this, we have gathered a comprehensive set of quasars at $1.5<z<7.5$ based on various samples of optical and near-infrared spectroscopic observations. Crucially, we take advantage of QUICFit (Meyer et al., 2019), a lightly user-supervised, model-independent spline fitting method for the quasar continuum and broad rest-frame UV lines redwards of Lyman-$\alpha$. This approach does not attempt to recover the systemic redshift, rather it simply retrieves in a uniform manner the observed peak of the BELs without assuming a particular model for these emissions. We can thus study the relative shift between C IV  and Mg II, but also the relative shifts between O I, C II, C III], Si IV  and the two previous lines as other potential markers of evolution. In doing so, we can probe the shifts of high-ionisation lines with respect to low-ionisation lines across cosmic time to evaluate on possible quasar evolution and provide insight into BLR physics.

The outline of this paper is as follows. In Section 2, we first detail our different datasets and the selection procedure of luminosity-matched control sample across our redshift range. We then assess the quality of our fits, and the different biases that could affect our measurement. Finally, we compare our method to PCA and Gaussian fitting on SDSS quasars. We present in Section 3 the evolution of the relative velocity shifts of all aforementioned BELs with respect to one another. In Section 4, we discuss our results in light of previous studies and we propose a novel interpretation of the increased C IV  blueshift observed in early quasars.

Throughout this paper, we adopt a $\Lambda$CDM cosmology with $\Omega_{M}=0.3$, $\Omega_{\Lambda}=0.7$, $H_{0}=70\text{ km s}^{-1}\text{Mpc}^{-1}$. We adopt a positive sign for velocity redshifts and negative for blueshifts.

2 Methods

2.1 Observational samples

In order to study the velocity shifts of quasar UV BELs across cosmic time, we have gathered a large sample of quasar spectra in infrared and optical light from different surveys and archival data. We briefly describe here the original samples from which our quasars were drawn and how we construct our luminosity-matched control samples. The relevant parameters of the studied datasets are summarized in Table 1 and the luminosity distributions are presented on Figure 1. We also show the stacked spectra of all our samples on Figure 2. Throughout this study, we refer to the specific luminosity at $1350$ Å, which we compute from the specific flux at $1350$ Å  rest-frame and the luminosity distance, simply as the luminosity. We note that although the $1350$ Å  luminosity is a practical choice given it can be measured for nearly all of our quasars, a different choice of luminosity matching (e.g. rest-frame X-ray or IR) might affect the results of this study due to dust correction effects. However, the emissivity at $1350$ Å  is more closely related to the UV emission lines of interest. In addition, the overlapping spectral coverage for our large sample of quasars is limited to the rest-frame UV.

SDSS DR7 and DR12: SDSS-I/II, -III (York et al., 2000; Eisenstein et al., 2011) delivered nearly $5$ million spectra and $32000$ deg${}^{{}^{2}}$ of imaging in the ugriz filters. In particular, we are interested here in two major catalogues of quasars observed and identified within SDSS data: the DR7 quasar catalogue (DR7Q, Schneider et al., 2010) and the DR12 quasar catalogue (Pâris et al., 2017). For DR7Q, we make use of the value-added catalogue of Shen et al. (2011) which provides additional properties for the quasars listed in the original DR7 catalogue. Shen et al. (2011) report redshifts derived from C IV, C III]  and Mg II, and the specific luminosities. In DR7Q, the line-based redshifts (i.e. redshifts derived from one emission line) are derived by fitting a multi-Gaussian template to each broad emission line and taking the template peak to derive the redshift. In DR12Q however, all quasars were inspected by eye and have thus a visual redshift, which is different from the line-based redshifts. The line-based redshifts in DR12 are derived from the peak of a PCA template fitted locally for each emission line (Pâris et al., 2011). There is some overlap between DR7Q and DR12Q, although the latter sample is larger and contains quasars with slightly lower luminosities and in a broader redshift range. The SDSS optical spectrograph (Smee et al., 2012) covers the range $3600-1000$ Å, with a resolution going from $R\sim 1500$ to $R\sim 2500$. Hence the detection of both C IV  and Mg II  is limited to quasars at $1.5\lesssim z\lesssim 2.4$. However, other UV BELs like O I, C II  or Si IV  can be detected simultaneously with C IV  in the complementary and overlapping redshift range $1.8\lesssim z\lesssim 5.1$.

XQ100: XQ100 (López et al., 2016) is a legacy European Southern Observatory (ESO) survey which observed 100 quasars at $3.5\leq z\leq 4.5$ at medium resolution ($R\sim 6000-9000$) and a high median Signal-to-Noise Ratio $SNR\simeq 30$ with the XShooter instrument (Vernet et al., 2011) at the Very Large Telescope (VLT). Its primary aims included absorption lines studies, AGN science, Broad Absorption Lines (BAL) studies, Damped Lyman-$\alpha$ absorbers (DLAs) and Lyman-$\alpha$ forest-based cosmology. The combined wavelength range of XShooter’s UVB, VIS and NIR arms ($3150-24800$ Å, or $3150-18000$ Å, depending on at which point in the survey the object was observed) is well suited to the aim of this study as all UV BELs can be retrieved over the entire redshift range $3.5\leq z\leq 4.5$, except when they fall in telluric regions, spectrograph arm overlaps and parts impeded by the atmospheric transmission. XQ100 also overlaps partially with SDSS, making this subset ideal to study the impact of resolution on the observed spectra and the recovered fitted continuum. However, the XQ100 quasars were selected to be intrinsically bright, making them difficult to luminosity-match with some of the high-redshift samples.

Giant Gemini GMOS Survey: The Giant Gemini GMOS survey (GGG, Worseck et al., 2014) obtained spectra of quasars over $4.4\leq z\leq 5.5$ with the Gemini Multi Object Spectrometers (GMOS) to study the Lyman-continuum (LyC) flux at $z\sim 5$ and determine the LyC photon mean free path. The publicly available data release of the GGG survey is composed of $163$ quasars at low-resolution ($R\sim 800-900$) with a medium-high $\text{SNR}\sim 20$, making it also suitable to capture the broad shapes of the UV BELs. The GMOS instrument covers the range $4800\AA\lesssim\lambda\lesssim 10200\AA$, meaning that C III]  and Mg II  can only be observed for quasars at $z\leq 4.4$ and $z\leq 2.6$, respectively. Thus we cannot measure the C IV-Mg II  relative shifts. However, since we do not restrict ourselves to this particular pair of BELs, we can use GGG quasars to compute velocity shifts of lower rest-frame wavelength BELs. GGG thus suitably completes the redshift ladder we have built by filling the gap between the SDSS datasets, XQ100 quasars and the highest-redshift samples at $z>5.5$.

$z\sim 6$** sample: ** Our first high-redshift sample (hereafter ‘z6’) consists of a subset of quasars from the extensive list of $z\sim 6$ object from Bosman et al. (2018), of which some where re-reduced later by Meyer et al. (2019). We keep only those with near-infrared coverage and high SNR, limiting the sample to $11$ objects. This approach ensures a comparison as complete as possible over $1.5\leq z\leq 6.5$, including all the lines from O I  to Mg II. This sample is critical to probe with additional spectra claims in the literature of (non-)evolution from the C IV-Mg II  velocity shifts. The $11$ quasars selected have both the highest-resolution ($R\sim 9000-11000$) and highest median SNR ($90$) of our samples (see Table 1).

$z\sim 7$** sample: ** Our second high-redshift sample is composed of the Mazzucchelli et al. (2017) quasars that are observed far enough in the infrared to detect at least C IV, to which we add the present record-holding $z=7.54$ quasar (Bañados et al., 2018), kindly provided by private communication of the authors. This sample thus comprises $12$ quasars observed with various instruments with $\langle z\rangle=6.73$, median $\text{SNR}\sim 6$ (although individual spectra range from $\sim 2$ to $\sim 100$) and resolution $R\sim 6000$.

Since BEL shifts correlate with luminosity, studying their potential evolution with redshift requires doing so at a fixed luminosity. We therefore construct luminosity-matched control samples from the SDSS quasars for all our higher-redshifts samples by oversampling their luminosity distribution. For for each quasar in a test high-redshift sample (XQ100, GGG, z6, z7), we draw a random SDSS quasar with a similar luminosity ($\Delta L=0.2$ dex). We then repeat the procedure until we reach $5000$ SDSS quasars. DR7Q and DR12Q only contain information about the peak location of C III], C IV  and Mg II, derived with different methods highlighted above which cannot be compared directly with our measurements. To enable a fair and uniform comparison, we also construct a small luminosity-matched (to z6) control sample of $108$ quasars spanning the redshift range $1.5<z<4.5$. This smaller control sample can be then fitted with our spline technique to control our systematics between different peak retrieval methods, resolution and instruments, and derive O I, C II  and Si IV-based redshift for SDSS quasars as well.

2.2 Continuum fitting and line-derived redshifts

We now describe the continuum fitting method and the subsequent derivation of the line-based redshifts that we apply uniformly to the quasar samples described above.

To enable the algorithm to run on a globally flat spectrum, each quasar is fitted with a continuum powerlaw (e.g. Vanden Berk et al., 2001) with the following form

[TABLE]

where $F_{0}$ [$10^{-16}$ erg s*-1* cm*-2* Hz*-1*] is the specific flux at $2500$Å, $\lambda$ the rest-frame wavelength in Å  and $\alpha$ is the powerlaw exponent. The powerlaw is fitted to regions relatively devoid of features following Decarli et al. (2010); Mazzucchelli et al. (2017) at $1285-1295,1315-1325,1340-1375,1425-1470,1680-1710,1975-2050,2150-2250$ and $2950-2990$ Å  rest-frame with a positive linear prior on $F_{0}$ and $\alpha$.

Once each rest-frame UV spectrum spectrum is divided by the powerlaw, we fit it using QUICFit, which we describe briefly here. Further details can be found in Meyer et al. (2019). QUICFit is a lightly-supervised spline quasar continuum algorithm which aims to avoid the human selection of appropriate continuum portions of the quasar spectra. Instead, the algorithm first discards BALs (typically of equivalent width EW$>1-2$ Å) by running a Gaussian matched-filter with a width of a few Ångströms. Smaller emission and metal absorption features, as well as bad pixels and cosmic rays are removed by searching for portions of the spectrum with excessive variance in the pixel-to-pixel flux differences. Indeed, if the spectral resolution is much higher than the continuum typical variation scale, the pixel-to-pixel flux increments and decrements distribution should follow the error array distribution. By computing the empirical variance in small portions of the spectrum, QUICFit can efficiently detect and removes structures (such as narrow absorbers, emitters or cosmic rays) that differ from a smooth continuum and expected uncertainty. Once this is done, the remaining parts of the spectrum can be considered representative of the intrinsic continuum and used to fit third order splines. We subsequently call these selected parts of the spectrum ‘continuum pixels’ later in the paper.

The spline knot points (where the polynomials of the spline join and the continuity conditions are enforced) can be initialized by the user, but it is found that one knot point at the location of each BEL, and one intermediate point in between, are sufficient to fit adequately third-order splines on the pre-selected continuum pixels. The process is lightly supervised and fast (a minute or two per quasar), as the user only has to adjust a threshold to reject pixels in narrow absorptions based on the expected uncertainty (see above) and choose whether or not to fit faint BELs by adding more spline knots. It is also found that complex lines such as C IV  or C III]  are better fitted with an extra knot point. Therefore we take particular care in accurately fitting the emission lines of interest: O I, C II, Si IV, C III], C IV  and Mg II. This approach avoids any misfits occurring in pipelines for large datasets such as SDSS quasars, while being still applicable to smaller datasets with a range of SNR and resolution.

The resulting spline fits and the line-based redshift solutions for all the BELs of our quasars are presented in Appendix A. Before investigating the suitability of our scheme to recover unbiased measurements of the BELs peak, we first focus on the quality of the fits by looking at the residuals. In order to do so, we assume that our selection of continuum pixels is correct, and compute the stacked residuals on these points only. We present on Fig. 3 (black points) the residuals of the continuum fits for all the different samples. The residuals are flat and show no particular trend, evidencing that we fit correctly the continuum of the selected continuum pixels.

Once the BELs have been fitted, we retrieve the peak of the spline fit to derive the line-based redshifts. In doing so, we take care to select only significant relative maxim in the continuum according to the following criteria. We first take the maximum in a rest-frame window of $10$ Å (O I, C II, Si IV, C III]) or $40$ Å  (C IV, Mg II) around the expected BEL location, computed using the reported redshift of the quasar in the literature. We then require that this maximum is significantly higher than the surrounding continuum. In order to do so, we take the relative minima of the continuum at lower and higher wavelength and we interpolate the continuum between the two minima. We compute the equivalent width (EW) of the line with respect to the interpolated continuum and we retain the line if it has $\text{EW}>3$Å. As the EW of the C IV  line anti-correlates with its blueshift (e.g. Coatman et al., 2016), we expect such a threshold could potentially bias the measurement of the blueshifts towards lower values. However this should not affect a potential evolution of the average blueshift. If multiple such maxima are retained, as can be the case for the Si IV  and C III]  blended multi-line complexes, we retain only the lower and higher wavelength peak, respectively. This is because we are interested in the Si IV  and C III]  emissions, not the blended O IV $\lambda 1402$ Å  and Si III] $\lambda 1892$ Å, respectively. This procedure disentangles the peaks when the two emission lines have similar equivalent widths and are resolved. Nonetheless, in some cases the line-based redshift simply tracks the peak of the blend, and is thus slightly biased towards larger (lower) values for Si IV  (C III]). We further comment on the impact on our results and the suitability of these line complexes to study relative velocity shifts in Section 3.

2.3 Instrument and resolution biases

.

Our quasar samples (Table 1) have different resolutions and SNR, and have been observed at different redshifts. We have fitted the continuum with QUICFit in a uniform way for the SDSS, XQ100, GGG, z6 and z7 samples. However, we need to ensure that our method is not biased by the varying resolution of the different samples towards lower or higher velocity shifts. This might be the case in the presence of a large number of unresolved associated narrow absorbers that would blend with the intrinsic profile and damp the blue or the red wing of the line. At higher resolution, these absorbers are resolved and thus successfully removed from the continuum to recover the intrinsic profile. This might create a marked difference between the high- and low- resolution samples if absorbers are found preferentially in the blue or red wing of BELs.

In order to address this point, we have degraded the XShooter spectra to SDSS-like resolution by re-binning them by a factor 3-5 (depending on the XShooter arm and the redshift of the object) and re-fitted the degraded quasar spectra with QUICFit. We demonstrate qualitatively the resilience of the fitting method by comparing the BEL spline fits recovered from the high- and low-resolution spectra in Figures 3 and 4. The recovered continua and residuals for O I, C II,Si IV  and C III]  are essentially identical. Nonetheless, the recovered continua for C IV  and Mg II  is slightly lower at the peak of the BEL. As expected, the lower resolution smoothes out sharp features such as the line peaks and narrow absorbers. This causes the residuals of the low-resolution continuum to be on average $1.0\sigma-1.5\sigma$ lower at the peak value. However the peak dampening of C IV  is highly symmetric, while Mg II  presents a small offset to redder values. We now quantitatively qualify the potential small bias induced in the final measurement of the line-based redshifts and subsequent relative velocity shifts.

We compare the redshifts derived from the fitted BELs at low and high resolution and we find that the line-based redshifts are not biased when the resolution is lowered (Fig. 5 and 6). The average difference between high- and low- resolution is smaller than natural intrinsic scatter of the measured redshifts by a factor of $5-10$ (see Table 2). Specifically, the mean redshift error due to the resolution decrease is $<100\text{ km s }^{-1}$, whereas the standard deviation is of the order of $\sim 400\text{ km s }^{-1}$.

To ensure that a change in both resolution and instrumental setup does not impact the measurement of the line-based redshift and relative velocity shifts, we apply our method to the $53$ XQ100 quasars observed with the SDSS spectrograph. We fit the quasar continuum with QUICFit directly on the SDSS spectra and retrieve all possible BELs when they fall in the spectrograph coverage. We present as before the scatter in the measurement of the line-based redshifts on Figure 7 and the subsequent impact on the derived velocity shifts on Figure 8. The values of the mean, median and standard deviation of the line-based redshifts errors are presented in Table 2 alongside the previous results comparing the original XQ100 XShooter spectra to the degraded resolution version. As expected, we find no further bias by comparing XShooter and SDSS spectra than by degrading XShooter spectra. The standard deviation of the velocity shift difference between the XShooter and SDSS is often smaller than for the previous case, which we attribute to being pessimistic whilst degrading the resolution. We find no bias either in the line-based redshifts or the derived BELs relative velocity shifts. We are thus confident that our spline fitting and peak retrieval method will not bias the measured relative velocity shifts between broad lines and is independent of resolution, instrument and observational setup and date.

2.4 Comparing QUICFit with PCA and Gaussian template fitting methods

QUICFit provides a resolution- and instrument- resilient method to fit the UV continuum and BELs of quasar continua. In this section, we compare line-based redshifts derived with our method to PCA and Gaussian fitting techniques used in previous studies. In doing so, we can also check that a change in fitting method does not bias the measurement by examining XQ100 quasars also observed in SDSS. SDSS DR12 is too large a sample to be fitted by hand, and the DR12 quasar catalogue only contains redshifts for C IV, C III]  and Mg II  (Pâris et al., 2017). These redshifts are derived from the peak of the PCA template (Pâris et al., 2012). DR7 redshifts, however, are derived using a Gaussian template fitting method (Shen et al., 2011). Both methods have their drawbacks and advantages, but they agree on the average C IV-Mg II  velocity shift at $1\lesssim z\lesssim 3$. In order to harness the statistical power of the DR7Q/DR12Q catalogues which measured C IV-Mg II  velocity shifts, we have to assess the possible difference between our measurements and the SDSS DR7/DR12 line-based redshifts. In order to do so, we compare the velocity shift derived from our method to the values provided by DR7Q and DR12Q where applicable in our luminosity-matched SDSS sample (Fig. 9). We show that there is no significant bias between our values, DR7Q and DR12Q values. The standard deviation between the different methods is however quite large. We find that the velocity shift difference between our results and DR7Q is ($\Delta v=(-83\pm 635)\text{ km s}^{-1}$ and $\Delta v=(47\pm 475)\text{ km s}^{-1}$ ($\sim 500\text{ km s}^{-1}$) for DR12Q. In both cases, the standard deviation is greater or consistent with the typical error expected from a change in resolution or a different instrumental setup as found in Section 2.3.

Finally, we can check that a concurrent change in resolution, SNR, instrument, and fitting method does not affect the measurement. For quasars present in both XQ100 and DR12Q, we compare the C IV- and C III]-derived redshifts from DR12Q to redshifts computed from our spline fits to the XShooter spectra. Because XQ100 contains quasars at $3.5<z<4.5$, the corresponding SDSS observations do not cover Mg II. We present on Fig. 10 the distribution of the two methods’ redshifts differences. We find that our solutions are similar to the SDSS ones with the usual $\sim 500$ km s*-1* spread for the two BELs species. The median shift, mean shift and the std deviation of the velocity shift distribution are $\overline{\Delta v_{\text{C{\tiny IV}}}}=53$ km s -1, $\langle\Delta v_{\text{C{\tiny IV}}}\rangle=90$ km s -1, $\sigma_{\Delta v_{\text{C{\tiny IV}}}}=392$ km s -1, and $\overline{\Delta v_{\text{C{\tiny III}}]}}=-7$ km s -1, $\langle\Delta v_{\text{C{\tiny III}}]}\rangle=-178.6$ km s -1, $\sigma_{\Delta v_{\text{C{\tiny III}}]}}=652$ km s -1 respectively. There are $\sim 10$ % of spectra presenting shifts above $500$ km s*-1* for C IV  and $\sim 30$ % for C III].

For these ‘catastrophic’ outliers, we plot their SDSS spectra, the XShooter spectra and the DR12Q and QUICFit line-based redshift solutions on Figures 16 and 17. Upon further inspection, we find that in most cases PCA seems to fail when multiple narrow absorptions, resolved with XShooter but not SDSS, artificially damp one wing of the emission. One could argue that PCA actually recovers the original feature because it is trained on unabsorbed spectra. However, the PCA emission line redshift solutions are determined by fitting a few eigenspectra to each UV broad line individually (Pâris et al., 2012). The PCA emission line redshifts are thus not derived from the overall spectral PCA fit and thus cannot reconstruct the ‘true’ emission from other features but rather perform a local best fit to the the emission shape. Another feature that the PCA seems unable to capture is a prominent contribution of Si III] $\lambda 1892$ Å  to the C III]  complex.

3 Relative velocity shifts of broad UV emission lines

We have demonstrated that our method is suitable for application to a large compilation of quasars of different SNR, resolution and redshift. We have also shown that our method is comparable to two main methods for deriving line redshifts (PCA templates or Gaussian fitting), with an intrinsic spread between methods of $\Delta v\sim 500\text{ km s}^{-1}$. We can now harness the statistical power from our extended sample as well as our extended number of studied BELs.

3.1 Single species broad lines: O I, C II, C IV  and Mg II

We start with the relative velocity shift of C IV  with respect to Mg II, as it is the focus of most previous studies. We note that only a fraction of our z6 and z7 quasars can be used for this analysis, as it requires good coverage of both the Mg II  and C IV  regions. The GMOS wavelength coverage does not allow us to detect the Mg II  line, therefore we discard the GGG sample for this particular analysis. We present on Figure 11 the results for both the z6 and z7 samples compared to appropriate DR12Q control samples. We find that, when luminosity is matched, there is a significant evolution of the shifts at $z\gtrsim 6-6.5$. Specifically, the mean C IV  blueshift evolves from $\sim 1000\text{ km s}^{-1}$ at $z\lesssim 6$ to $\sim 2500\text{ km s}^{-1}$ at $z\sim 7$, as already reported by Mazzucchelli et al. (2017). The double-sided KS-test, which tests whether two samples are drawn from the same different distribution, gives no evidence ($p=0.69$) for the z6 blueshifts being different from their luminosity-matched DR12Q samples. This result is similar to that reported by Shen et al. (2019) . However, the p-value for the z7 sample ($p=3\cdot 10^{-6}$) is clear evidence for evolution compared to luminosity-matched SDSS objects. The evolution of the C IV  line is also clearly visible in the stacked spectra of all our samples (Fig. 2), despite the intra-object spread strongly smoothing the feature in the stacked spectrum.

We now turn to the relative velocity line shifts for four BELs (O I, C II, C IV, Mg II) in all our different samples. In order to do so, we compute the mean relative shift for every pair of species in all samples, which we show on Figure 12. The errors are computed by boostrapping using $1000$ samples of the size of the smallest sample (z6, $N=11$) to make them comparable across redshift and samples. The low-ionisation lines display no significant velocity offsets either between themselves at any redshift. The small relative blue- or red-shift of the lines are all consistent with bootstrapping errors. These low-ionisation line are not predicted to shift with respect to each other, therefore an absence of relative shifts is a good indicator of an unbiased measurement of the line shifts.

The C IV  line, however, is markedly blueshifted from all the lower-ionisation lines (O I, C II, Mg II) at all redshift. The mean blueshift significantly increases at $z\gtrsim 6.5$ by a factor $2.5$, irrespective of the choice of low-ionisation species used as a reference to measure the shift of C IV. This change is indicative of different environments for the emission (and possibly the absorption) of C IV  and low-ionisation lines. It is probably linked to a dichotomy between high- and low- ionisation lines given that they are thought to originate in different parts of the BLR. To test this hypothesis, measuring the mean blueshift of another high-ionisation line is needed. However, the only other high-ionisation line in the rest-frame UV is N V $\lambda 1240$ Å, which is often blended with Lyman-$\alpha$ and low-ionisation Si II $\lambda 1260$ Å. Moreover Lyman-$\alpha$ is increasingly absorbed at high-redshift, making the deblending more difficult. Another potential high-ionisation line, He  II $\lambda 1640$ Å  is often too faint to be detected in most high-redshift quasars. Hence, spectroscopy of the rest-frame optical will be needed to conclude whether this blueshift is a particularity of the C IV  line or a dichotomy between low- and high-ionisation lines.

3.2 Broad lines complexes: Si IV  and C III]

We have left aside the Si IV  and C III]  broad lines from the previous analysis. These two intermediate-ionisation lines are known to be complexes: Si IV $\lambda 1397$ Å  is blended with O IV $\lambda 1402$ Å  whereas C III]  $\lambda 1909$ Å  is blended with Al III]  $\lambda 1857$ Å  and Si III]  $\lambda 1892$ Å. These complexes could affect the measurement and the interpretation of any potential velocity shifts, especially since the relative line strengths of the blended emission lines appear – qualitatively – to be evolving with redshift. Figure 2 illustrates this effect: the O IV line disappears between $z\sim 5$ and $z\sim 6$, whereas the Si IV  BEL increases in strength at high-redshift, even though the samples are luminosity-matched. Notwithstanding this complication, we attempted to measure the velocity shifts of these complexes by taking the red/blue peak of the complex accordingly when the feature was double-peaked, and the overall peak of the complex when the lines were blended (see Section 2).

We find no evidence of velocity shifts of C III]  and Si IV  relative to low-ionisation lines, and the C IV  line is seen to blueshift with respect to these lines in the same manner as to single low-ionisation lines (see Fig. 18). This result can be interpreted in two ways. If we were successful at disentangling the complexes, our findings indicate that C III]  and Si IV  are not blueshifted with respect to low-ionisation lines at any redshift. Even if we were not able to satisfactorily separate the blending, our results still provide some constraints on the velocity shifts of C III]  and Si IV. As C III] is dominant over Al III]  and Si III] at late times and those blended lines have shorter wavelengths, a relative strengthening of Al III]  and/or Si III] could only lead the peak of the blend to blueshift. Redshifting of the complex peak would then be attributable to C III]. Since neither is seen, we can plausibly conclude that on average C III] does not redshift with respect to low-ionisation single species. Moreover, any intrinsic blueshift would need to be carefully masked by evolution of the other blended lines in order to go unnoticed. The situation is reversed for Si IV, which is blended with lines of longer wavelength than itself. The lack of velocity shift strongly suggests no blueshifting of Si IV, with the only alternative being a simultaneous equal shifting of O IV coupled with a reversal in line strength balance.

Previous studies have already pointed at the lack of a strong Baldwin effect and blueshifts of Si IV  (Osmer et al., 1994; Richards et al., 2002). If the similar lack of blueshift at high-$z$ is confirmed, this would imply that C IV  is seen to blueshift with respect to both low- and intermediate-ionisation lines, suggesting a different origin of the line within the BLR.

4 Discussion

4.1 Comparison to previous works

The most striking result of this analysis is the rapid evolution in the velocity blueshift of C IV  in the restricted redshift range $5<z<7$ depicted in Figure 12. Possible evolutionary trends in velocity shifts, specifically in the blueshift of C IV  with respect to Mg II, have been previously examined by several authors. Shen et al. (2019) found no significant evolution to $z\sim 6$ consistent with our results over $1.5\lesssim z\lesssim 6$. If we restrict an evolutionary test to luminosity-matched samples drawn from the z6 and SDSS datasets, the probability that these are drawn from two different distributions has an inconclusive p-value of $0.69$ (see Fig. 11). Although the mean velocity offset at $z\simeq 6$ is slightly larger ($\Delta v_{\text{C{\tiny IV}}-\text{Mg{\tiny II}}}=(-1400\pm 334)\text{ km s}^{-1}$), given the uncertainty due to the small sample size there is no convincing case for evolution.

However, with the addition of the $z\sim 7$ quasars from Mazzucchelli et al. (2017) and Bañados et al. (2018), we recover the larger mean blueshift of C IV  found by these authors. Since we have applied the same method to different samples over the redshift range $1.5\leq z\leq 7.5$, we argue this evolution cannot be caused by resolution or other instrumental effects. Our method also allows us to extend the analysis by substituting O I  or C II  for Mg II . Given the similarity of the trends so revealed, we can aggregate over the low-ionisation lines to construct a C IV-low-ionisation velocity shift versus redshift trend that uses all datasets (Fig. 13). When aggregating the low-ionisation species, we first construct a common ‘systemic’ redshift from the mean of all available low-ionisation line-based redshifts. We then derive the C IV  blueshift and compute errors by bootstrapping as before. The sharp increase in the C IV  blueshift between $z\sim 5-7$ is evident. We also show that the recently published C IV-Mg II  velocity shifts for individual $z\gtrsim 6.5$ quasars are in agreement with the increasing mean C IV  blueshift trend (Wang et al., 2018a; Pons et al., 2019; Reed et al., 2019).

There has been some debate on the possibility that biases in the line-based redshift determination of SDSS quasars might remove the most blueshifted C IV  lines. DR7Q derives C IV  blueshifts up to a maximal value of $5000\text{ km s}^{-1}$, making the comparison with some $z\sim 7$ quasars impossible. Coatman et al. (2016) showed that using a definition of the line center based on the peak emission or line centroid was strongly affecting the SDSS distribution of C IV  blueshifts, potentially alleviating some concerns raised with the first reports of large blueshifts in early quasars (De Rosa et al., 2014). The recent results of Reed et al. (2019) showcases the use of a currently unpublished Independent Component Analysis (ICA) determination of the SDSS C IV  blueshifts that identifies multiple extreme blueshifts $>4000\text{ km s}^{-1}$ in SDSS. Despite both studies extending the range of existing SDSS C IV  blueshifts, the mean shift and the abundance of redshifted C IV  lines still does not match that for $z\sim 7$ quasars. We have tested this by independently measuring the line-base redshifts of C IV  and different low-ionisation BELs in about a hundred luminosity-matched SDSS quasars. This result suggests that improvements in quasar template fitting will most likely not result in a dramatic increase of the average C IV  blueshift at low redshift. Nonetheless, these advances will enable the selection and study of low redshift quasars with extreme C IV  blueshifts, which we argue below will be a crucial point to further our understanding of reionisation quasars.

4.2 Interpreting the high-redshift increased mean C IV  blueshift

The most intriguing aspect of the trend in Fig. 13 is that the inferred evolution occurs in less than 1 Gyr. At low redshift, discussion of the nature of C IV blueshifts has centered around radiation-driven winds originating in the central regions, or the BLR, of SMBH accretion discs. The broadening of quasar UV emission lines was first modelled through primarily planar outflows, in which material is stripped from the accretion disc and accelerated outwards by X-ray radiation and/or resonant line radiation pressure from the inner regions (Krolik & Begelman, 1986; Murray & Chiang, 1995, e.g.). Based on the singular behaviour of the C IV BEL, Wills et al. (1993) first suggested a possible separate origin of the line in a hotter, polar wind component (see also Denney 2012). Later models showed how gas could be entrained at a steeper angle from more central regions (e.g. via increased magnetisation of the outflows, Proga et al. 2000; Elvis 2000) and even give rise to outflows which are significantly polar by entraining gas from a dusty torus, rather than a thin accretion disc (e.g. Gallagher et al. 2015). Early results from sub-millimeter imaging have supported the possibility of significantly perpendicular winds, with detections of elongated dusty emission around AGN in Seyfert galaxies being preferentially parallel to their ionisation cones, rather than perpendicular (López-Gonzaga & Jaffe 2016 and references therein). If low-ionisation and high-ionisation broad lines do indeed arise through such related but different processes along different axis, with higher-ionisation outflows being more polar, it is possible to conceive scenarios in which one evolves but not the other.

The most commonly invoked explanation for the fast evolution of C IV blueshifts at early times involves (unusually) strong BLR outflows or winds to explain (extremely) blueshifted C IV  emission (e.g. Richards et al., 2011; De Rosa et al., 2014; Mazzucchelli et al., 2017). This claim is supported by the widespread presence at $z\sim 7$ of extreme blueshifts with $\Delta v\gtrsim 3000$ km s -1 which are exceedingly rare at $z<5$. However, it is unclear what mechanism could power such an increase in wind speed in early quasars, or whether this can account for the whole effect. Such extreme blueshifts are found to be quite common in some lower-redshift populations such as the WISE/SDSS selected hyper-luminous quasars (WISSH, Bischetti et al., 2017; Vietri et al., 2018) at $z\sim 2-4$, and hence are not exclusive to $z\gtrsim 7$ quasars. Specifically, WISSH quasars with weaker O III $\lambda 5007$ Å  exhibit larger blueshifts and lower X-ray to optical luminosity ratio, which are also found in some high-redshift quasars (Bañados et al., 2018), and this trend seems to persist even in moderately bright SDSS quasars (Coatman et al., 2019). This suggests that mechanisms which increase C IV blueshift are still in place at much later times, and even dominant in some populations. Moreover, anomalously large C IV blueshifts at early times do not seem to correlate with higher accretion rates – as one might perhaps expect if they were caused primarily by anomalously strong winds. Mazzucchelli et al. (2017) reported that their $z\gtrsim 6.5$ sample had high accretion rates near the Eddington limit and their SDSS luminosity-matched sample did as well. Finally, a wind-only model would imply that some cases of no blueshift or even redshifted lines should always occur due to orientation with respect to the observer, unless the wind is increasingly spherical. In contrast, there is a total absence of redshifted C IV  emission lines in luminous quasars beyond $z\gtrsim 6$. If increased wind speeds do power the increase in C IV blueshifts, they are then most likely accompanied by a change in morphology. This conclusion either suggests a fundamental change in the properties of quasars as they emerge from the reionisation era or perhaps some selection effect which would cause quasars with non-blueshifted C IV  lines to be absent at $z\gtrsim 6$ despite our attempts to create luminosity-matched samples.

4.2.1 Orientation selection bias

In a simple ‘opaque torus’ model of quasars (Denney, 2012, and references therein), C IV  radiation originates either in the broad line region or a surrounding intermediate line region and launches fast winds perpendicular to the plane of the accretion disk (or at least significantly more polar than low-ionisation winds). If the accretion disk is viewed side-on where the line of sight is opaque, this will result in a broad C IV  emission with no blueshift and a dimmed UV quasar luminosity. Viewed increasingly face-on, the C IV  line will be blueshifted with a shape dependent on the wind velocity profile. This geometric model explains many observational features of broad quasar emission lines, such as the consistent, but moderate average blueshifts at low redshift.

Assuming that high-redshift quasars are obscured to the same extent as those at lower redshifts, the enhanced C IV  velocity offsets would then likely be related to the relative youth of quasars at the end of the reionisation era. Although the $z>6.5$ quasars in our sample have already massive ($\sim 10^{9}M_{\odot}$) black holes accreting at nearly the Eddington rate (De Rosa et al., 2014; Mazzucchelli et al., 2017; Bañados et al., 2018), the UV luminosity of quasars would be expected to be attenuated when they are viewed side-on, which may lead to them not being detectable at our luminosity threshold at $z\gtrsim 6.5$. In this picture, quasars with masses and accretion rates large enough to be UV-bright despite being side-on might not have yet assembled as early in cosmic time. Sampling the top of the UV luminosity distribution at any cosmic time would then unwittingly bias observations towards face-on objects with large C IV  blueshifts.

If this selection bias were the only effect at work, quasar continuum reconstruction models which successfully capture the variety of low-redshift quasars should be expected to work similarly well at high-redshift. Indeed, the physical origin of the C IV  blueshift – outflows and orientation – would be intrinsically the same across cosmic time. One might expect a range of related observational consequences, for example:

• •

The most extremely blueshifted objects of a given luminosity – being face-on – should present roughly the same $\Delta v_{\text{C{\tiny IV}}}$ at all redshifts. In our study, the maximum blueshift in the $z\sim 7$ sample is $2000\text{ km s}^{-1}$ greater than our DR12Q luminosity-matched sample extremum (Figure 9). As our samples are relatively small, it is possible that larger datasets of high-redshift quasars and new fitting techniques of low-redshift quasars could still change the values of the most extreme blueshifts.

• •

The shape and strength of the highly blueshifted C IV  emission lines should be the same across redshift, since the evolution of the mean shift would be due only to the disappearance of quasars with non-blueshifted C IV. We present in Fig. 14 the mean C IV  line profiles of our quasars at different redshifts, stacked according to the blueshift of the C IV  emission compared to low-ionisation lines. It can be seen that the higher-redshift quasars show a tentatively attenuated profile and smaller equivalent widths when C IV  is not blueshifted, as reported in previous studies (Mazzucchelli et al., 2017; Shen et al., 2019).

• •

The effect should vanish for fainter high-redshift quasars. Perhaps the best test of this ‘selection bias’ would be to examine the C IV  velocity shift at early times in less luminous objects. Such a strategy is obviously very challenging. However, the discovery of a faint $z\sim 7$ quasar with a significant C IV  blueshift ($\sim 2500\text{ km s}^{-1}$) in the Subaru High-$z$ Exploration of Low-Luminosity Quasars survey (SHELLQS) provides marginal evidence against this interpretation (Matsuoka et al., 2019) while a recent small sample of six relatively faint quasars at $z\sim 6.3$ show mixed results including potentially the first C IV redshift with respect to Mg II at $z>6$ (Onoue et al., 2019). Clearly a larger sample of fainter, possibly lensed quasars above $z\gtrsim 6$ can test this hypothesis.

4.2.2 Increased obscuration

A competing or complementary alternative is that dust obscuration is more intense in early quasars. In this case, those quasars viewed side-on will again be (even more) strongly de-selected at high-redshift, increasing the average blueshift of the C IV  line. Treister et al. (2011) reached the conclusion that black hole accretion is mostly obscured in the Early Universe by looking at the X-ray emission at $z\sim 6-8$ attributed to early quasars. Trebitsch et al. (2019) showed that increased quasar obscuration at high-redshift could perhaps be due to chaotic feeding and instabilities in the less mature host galaxies. The same authors also showed that early quasars accreting at higher Eddington ratios are preferentially surrounded by dense, obscuring gas and are only being visible with $M_{1450}>-23$ about $4$% of the time. Given that the highest redshift quasars must have accreted at nearly Eddington rate for most of their life (e.g. Bañados et al., 2018), it is thus likely that they are more obscured than the lower-redshift ones. By selecting the most UV-bright objects at high-redshift we would then only find quasars which have both very high accretion rates and are preferentially orientated, with potentially very narrow visibility channels. We speculate that chaotic feeding of the central black hole would affect the BLR and thus the BELs. Models and templates extracted from low redshift quasars would thus be expected to fare fairly poorly on high-redshift objects.

This second scenario might also explain some other intriguing recent finds in reionisation era studies. The small near-zones of some high-redshift quasars (Eilers et al., 2017; Eilers et al., 2018b) could be due to intermittent obscuration of the quasar, blocking LyC photons up to $z\sim 6$, where the near-zone would start to grow steadily and C IV  blueshift normalize. Depending on the type of obscuration at work, highly obscured and thus UV-faint quasars might in fact already be known as faint AGNs hosted in high-redshift galaxies showing N V emission lines (e.g. Stark et al., 2017).

Finally, both interpretations have potential consequences for the contribution of quasars to hydrogen reionisation (e.g. Madau & Haardt, 2015). In a ‘geometrical selection effect’ scenario, only quasars which are ‘sufficiently face-on’ can be currently detected at high redshift. Indeed, $\sim 50\%$ of z7 quasars have C IV  blueshifts $\geq 2500$ km/s, while those make up $<15\%$ (at $2\sigma$) of the DR12Q sample. Since quasar ionising radiation output is calculated over all angles, this would imply that an average high-redshift quasar in our luminosity-matched sample contributes less to the IGM photon budget than its low-redshift counterpart. We deduce this because edge-on, non-C IV-blueshifted quasars on average are expected to have larger masses and larger accretion rates in an orientation-only model – and those are the ones missing at early times. By contrast, in a scenario where early quasars are more obscured than their low-redshift counterparts, it is likely that we are currently observationally missing a large fraction. We would expect to find an increased number of obscured early AGN – which are perhaps contributing to the bright end of the galaxy UV luminosity function (Ono et al., 2018) or to the AGN X-ray luminosity function, which should be less obscured (Giallongo et al., 2015). Using the rough value of $\sim 4\%$ UV-bright visibility fraction from Trebitsch et al. (2019), and if obscuration is mostly geometric, we could be underestimating the number density of early massive quasars by factors of $3-4$. This deficit could be somewhat balanced by less ionising emission escaping through smaller channels. In any case, the C IV  blueshift conundrum has far-reaching consequences on the role of quasars in reionisation.

Acknowledgements

We authors thank the anonymous referee for the thoughtful comments on the manuscript. RAM, SEIB and RSE acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 669253). We would like to thank E. Bañados and C. Mazzucchelli for useful discussions and kindly sharing their data with us. We thank M. Trebitsch for useful discussions.

Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme(s) 084A-0550, 084A-0574, 086A-0574, 087A-0890 and 088A-0897.

Based on observations made with ESO Telescopes at the La Silla or Paranal Observatories under programme ID(s) 189.A-0424(A), 189.A-0424(B)

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Appendix A Line-derived redshifts for all quasars

We present the continuum splines fitted around the broad emission lines for each quasar in this study. Panels from left to right show broad emission lines arranged by increasing rest-frame wavelength, i.e. : O I, C II, Si IV, C IV, C III], Mg II. For each quasar and line, we show the observed flux (black), the fitted continuum splines (blue line) and the location of the peak if it satisfies the criteria of Section 2 (red vertical liens). For XQ100 quasars, we also show the peak recovered from the low-resolution version of the spectra (green vertical lines) and the one derived from the SDSS spectra if the object is present in SDSS (blue dotted vertical lines). Each BEL and fitted continuum are rescaled to facilitate the inspection. We thus caution that the lines are all at a different scale except for O I  and C II  which are rescaled by the same amount as they are contiguous in rest-frame wavelength. For the XQ100 quasars, we also show the continuum splines fitted to the resolution-degraded flux, and the derived peak positions thereof. The quasar name and literature redshift are indicated on the left of all panels. We give only the fits for the first $6$ objects of the XQ100 sample to limit the length of this manuscript. The remaining plots are available as online supplementary material.

Eventually, we present a summary table of our measured quasar redshift based on each BEL in Table 3. The entire tables for each samples are available as supplementary material.

Appendix B Catastrophic redshift errors between our method and DR12Q tabulated values

Appendix C Relative velocity shifts including Si IV  and C III]

We report on Fig. 18 the relative velocity shifts of all lines, including the Si IV  and C III]  complexes. It must be noted that, because shifts are derived from the peak of the complex, the shift of these two complexes are mostly driven by the relative amplitude and thus abundance of the underlying emission lines evolving with redshift.