© ESO 2020

1. Introduction

Ross et al. (2015) discovered a population of extremely red quasars (ERQ) in Data Release 10 (DR10) of the Baryon Oscillation Sky Survey (BOSS; Ross et al. 2012; Dawson et al. 2013) in the Sloan Digital Sky Survey-III (SDSS-III; Eisenstein et al. 2011). These ERQ show very red colours (SDSS r band to WISE W4 band) similar to dust-obscured galaxies (DOG). The authors identified 65 ERQ that spanned a redshift range of 0.28 <  z <  4.36 with a bimodal distribution, with peaks at z ∼ 0.8 and z ∼ 2.5. Most objects are type 2 quasars (QSO2) or heavily reddened type 1 quasars (QSO1), but a subsample of 12 objects presents very peculiar emission line to continuum properties that defy standard explanations based on extinction and/or orientation. This subclass, which was later called core ERQ by Hamann et al. (2017, hereafter H17), is the topic of this paper.

Hamann and collaborators enlarged the core ERQ catalogue to 97 objects by defining a less stringent colour condition. These 97 objects have nearly uniform peculiar properties selected through i − W3 ≥ 4.6 (AB system) and a rest-frame equivalent width (REW) of CIVλ1550 (CIV hereafter), REW_CIV >  100 Å, at redshifts z ∼ 2.0–3.4.

Core-ERQ have very high bolometric luminosities (median log(L_bol(erg s⁻¹)) = 47.1 ± 0.3; see H17 and Ross et al. 2015 for a detailed characterisation and investigation of core ERQ properties). They show unexpectedly flat UV spectra given their red UV to mid-IR colours and large line REW: 50% of the core ERQ have REW_CIV >  150 Å compared to ≪1% in normal blue QSO1 matched in i or W3 magnitude. They show signs of strong absorption in the X-rays, with inferred column densities of N_H ≳ 10²³ cm⁻² (Goulding et al. 2018).

Core-ERQ have peculiar, wingless rest-frame UV emission line profiles. The values of the full width at half-maximum (FWHM) lie between those found for very luminous type 1 and type 2 active galactic nuclei (AGN) at similar z. Their median FWHM_CIV = 3050 ± 990 km s⁻¹ is narrower than the FWHM_CIV = 5836 ± 1576 km s⁻¹ for blue quasars that were matched in W3 (H17), but significantly broader than the FWHM_CIV <  2000 km s⁻¹ of QSO2 (Alexandroff et al. 2013) and narrow-line radio galaxies (NLRG, De Breuck et al. 2000) at similar z. High blueshifts in excess of 2500 km s⁻¹ in CIV and other high-ionisation UV lines compared to the HI Balmer and low-ionisation permitted lines in the UV are also reported, which have been interpreted in terms of outflows (H17). CIV blueshifts are common in normal QSO1 (e.g. Gaskell 1982; Sulentic et al. 2007; Runnoe et al. 2014; Vietri et al. 2018). They have been widely interpreted within the context of accretion disc wind models (e.g. Richards et al. 2011)

The kinematic properties of the forbidden [OIII]λλ4959,5007 lines ([OIII] hereafter) revealed by near-infrared spectroscopy are even more striking (Zakamska et al. 2016, hereafter Z16; Perrotta et al. 2019). Core-ERQ exhibit extremely broad and blueshifted [OIII] emission, with widths (W₉₀) ranging between 2053 and 7227 km s⁻¹, and maximum outflow speeds (V₉₈) up to 6702 km s⁻¹. According to the results reported by Zakamska and Perrotta, at least 3–5% of their bolometric luminosity is being converted into the kinetic power of the observed winds, which means that this has the potential of affecting the entire host galaxy.

Core-ERQs also show peculiar line ratios. For instance, Lyα is often very strongly absorbed, sometimes almost completely. This results in very low line intensities relative to other lines such as CIV or NVλ1240 (hereafter, NV). Core ERQs are also characterised by a high NV/CIV ratio (often > 1.5) and an intermediate [OIII]/Hβ ∼ 1–4 ratio (Perrotta et al. 2019)¹ between type 1 and type 2 luminous AGN.

Distinguishing the various explanations for these puzzling spectral features of core ERQs is of great interest in its own right. Moreover, these systems might also be highly relevant to studies of galaxy formation and evolution. They have been proposed to be near- or super-Eddington accreting obscured quasars, which are hosts of some of the most massive black holes (BH) at z ∼ 3 that can trigger strong galactic outflows that inhibit star formation in the early universe. They may represent an intermediate phase in which a heavily obscured quasar blows out the circumnuclear interstellar medium (ISM) with very energetic outflows before they become an optical quasar (Z16; H17; Goulding et al. 2018; Perrotta et al. 2019).

In the unified model of AGN (e.g. Antonucci 1993; Urry & Padovani 1995), the orientation with respect to the observer of a dusty, obscuring central structure (torus or other) located within the BH gravitational radius of influence can explain certain differences found between an obscured and an unobscured quasar. Such an obscuring structure would block the view along some lines of sight towards the accretion disc and the clouds within the broad-line region (BLR), so that these become partly or totally hidden. Independently of the specific properties of the blocking structure and the role of other factors that can influence the diversity of quasar properties, it is clear that orientation is key for explaining certain differences (e.g. Alonso Herrero et al. 2011; Ramos Almeida et al. 2011; Elitzur 2012; Mateos et al. 2016). We investigate in this paper the role played by orientation, extinction, and the extreme luminosities of the core ERQ on their observed emission lines and on the properties of their continuum spectra, together with the possible implications on the evolutive scenario. We adopt H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.73, and Ω_m = 0.27.

The sample we studied consists of 21 core ERQ. We focus on the 20 core ERQs studied by Perrotta et al. (2019; which include the 4 objects studied by Z16), for which the authors presented near-infrared spectroscopy in the Hβ–[OIII]λλ4959,5007 region. These lines, especially the [OIII] doublet, will provide essential information in our argumentation. We excluded core ERQ-like objects (see H17) so that we obtained a more homogeneous sample in terms of peculiar spectral energy distributions (SED) and emission line properties. Detailed information on the sample can be found in Perrotta et al. (2019). We also included the core ERQ SDSS J171420.38+414815.7 at z = 2.34 (J1714+4148 hereafter) that we observed with the 10.4 m Gran Telescopio Canarias (GTC) (see Sect. 2).

We measured UV line ratios for the sample of core ERQ using the BOSS optical spectra. As a guidance, we show in Table 1 the line ratios for the four core ERQ in Z16 and for J1714+4148. The UV ratios for 16 the remaining core ERQ are shown in Appendix A. Throughout the text, we refer to the emission lines as follows: NV for NVλ1240, CIV for CIVλ1550, HeII for HeIIλ1640, CIII] for CIII]λ1909, and SiV+OIV] for SiV+OIV]λ1400.

Lyα is often blended with NV (see Fig. 18 in H17). This effect is less severe in core ERQ than in QSO1 because Lyα is often heavily absorbed and the lines are narrower. These two lines were clearly isolated in 16 of the 21 objects. A careful evaluation of the possible contamination of NV by Lyα was necessary in 5 objects with severe blending (see, e.g. spectra of J0834+0159 and J2215−0056 in Fig. 1 of Z16).

To estimate this effect, we proceeded as follows. Maximum contamination might be expected for the highest possible Lyα flux. That is, for no Lyα absorption and maximum flux relative to NV. As a reference to estimate this, we used the measured CIV fluxes and the typical Lyα/CIV ratios observed in AGN.

A range of values of Lyα/CIV ∼ 1–20 is measured in luminous type 2 AGN (e.g. De Breuck et al. 2000; Villar Martín et al. 2007; Alexandroff et al. 2013). The lowest values are a consequence of prominent Lyα absorption. The ratio is typically Lyα/CIV ≲ 6 in QSO1 (Zheng et al. 1997; Lusso et al. 2015; Telfer et al. 2002).

For each object we created artificial Lyα spectral profiles for which Lyα/CIV = 20. The Lyα FWHM and velocity offset relative to NV were set to be equal to those of the strongest emission lines (CIV, SiIV). This “unabsorbed” Lyα profile would result in a maximum possible contamination of NV. The spectral profiles of the strongest emission lines in our sample were reasonably well fitted with a Gaussian. Consistently, in our analysis we assumed a Gaussian shape for the unabsorbed Lyα. Because the UV lines of each target typically show a range in their FWHM and in their redshifts z_i, we assumed the worst-case scenario and adopted as Lyα FWHM the broadest value observed among the different lines. It was found that uncertainties on the relative velocity (z_Lyα) with respect to NV did not have a significant effect on the results.

The uncertainties of the NV flux derived from this evaluation were ≲20% for four of the five objects under consideration and ∼30% for a fifth object. Because Lyα is heavily absorbed in general and the worst-case scenario was adopted when we estimated the Lyα contamination, we can reasonably conclude that its effect is low in general. The errors hereafter quoted for the line ratios involving NV include all the uncertainties described above. They can therefore be qualified as the maximum expected errors.

2. J1714+4148: GTC EMIR observations, data reduction, and analysis

We obtained HK-band spectroscopy of the core ERQ J1714+4148 at z = 2.34 with the 10.4 m GTC and the EMIR (Espectrógrafo Multiobjeto Infra-Rojo) instrument in long-slit mode (program GTC14-19A). EMIR is a near-infrared wide-field imager and medium-resolution multi-object spectrograph installed at the Naysmith-A focal station. It is equipped with a 2048 × 2048 Teledyne HAWAII-2 HgCdTe near-infrared optimised chip with a pixel size of 0.2″. The low-resolution HK grism covers a spectral range of ∼1.45–2.42 μm with a dispersion of 6.8 Å pixel⁻¹. The slit width used during the observations was 1.0″. The instrumental profile measured from the OH sky lines is FWHM_IP = 31.3 ± 4.0 Å. The total exposure time on source was 3840 s. A typical ABBA nodding pattern was applied. The seeing during the observations was FWHM ∼ 0.7″.

The spectra were reduced using several python routines customised by GTC staff for EMIR spectroscopic data. The sky background was first eliminated using consecutive A-B pairs. They were subsequently flat-fielded, calibrated in wavelength, and combined to obtain the final spectrum.

To correct for telluric absorption, we observed a telluric standard star with the same observing set-up as the science target, immediately after the J1714+4148 observations and at similar airmass. To apply the correction we used a version of Xtellcor (Vacca et al. 2003) specifically modified to account for the atmospheric conditions of the La Palma observatory (Ramos Almeida et al. 2009).

Relative flux calibration was applied using the spectrum of the star, which was obtained with the same narrow 1.0″ slit. The accuracy is ∼10%. It is not clear how different this is from the absolute flux calibration because no near-infrared magnitudes are available for the object. Based on previous experience, we estimate a maximum deviation of ≲30% between the relative and absolute flux calibrations. Because galactic extinction is very low (A_V = 0.06), no correction for this effect was applied.

We fitted Hβ and [OIII]λλ4959,5007 with single Gaussian profiles to obtain an approximate quantitative characterisation of the gas kinematics. The separation in λ and the flux ratio (3:1) for the [OIII] doublet were fixed to the theoretical values. Each line was parametrised with a central λ, observed FWHM, and amplitude (and thus, the flux). The two [OIII] lines were forced to have the same FWHM. More complex profiles are likely to be more realistic (Perrotta et al. 2019), but this approach is not possible with our data because of the low signal-to-noise ratio (S/N; [OIII]λ5007 is detected with S/N ∼ 6). All FWHM values in km s⁻¹ quoted below were corrected for instrumental broadening in quadrature.

We show in Fig. 1 (first panel) the Hβ and [OIII] spectrum and the best fits for different sets of constraints. We attempted different fits by applying a variety of kinematic constraints to determine useful ranges of the line FWHMs and the shift in velocity Δv_{Hβ − [OIII]}. Some examples are shown in Fig. 1. Each panel shows the fit with the smallest mathematical errors (minimum reduced chi-squared, χ²) for that specific set of constraints.

Fig. 1. GTC EMIR spectrum of J1714+4148 in the Hβ and [OIII] spectral region (black lines). The flux is in units of ×10−18 erg s−1 cm−2 Å−1. Different fits (red lines) were attempted by applying a variety of kinematic constraints to [OIII] and/or Hβ (fits 2–7) or by leaving the kinematic parameters free (fit 1). For each fit, the constraints on the observed FWHM in Å are shown. It is also specified whether Hβ was forced to have the same redshift as [OIII]λ5007 (λ forced, fits 2 and 2b). The reduced chi-squared χ2 are also shown, and the best-fit profiles are presented underneath each panel.

The low S/N of the spectrum prevents accurate constraints on the [OIII] and Hβ parameters. The visual comparison between the data and the fits shows that χ²∼1 is not always associated with satisfactory results (see, e.g. fits 4 and 7).

Based on this analysis, we propose two tentative results: the [OIII] lines are broad, with a most probable range 1700 ≲ FWHM_[OIII] <  2300 km s⁻¹ ( Å. The fits suggest that Hβ is broader than [OIII] and is redshifted, although we cannot discard that this is an artificial result due to the noise of the Hβ profile. Tentatively, 2000 ≲ FWHM_Hβ ≲ 3000 km s⁻¹ and 900 ≲ Δv_{Hβ − [OIII]} ≲ 1300 km s⁻¹. Such broad Hβ and [OIII] lines and prominent [OIII] blueshifts are common in core ERQ (Perrotta et al. 2019).

3. Results and discussion

3.1. Comparison of the UV emission line ratios with blue QSO1, QSO2, and high z NLRG

The line ratios of the 4 core ERQ in Z16 and J1714+4148 are shown in Table 1 (see Appendix A for the remaining 16 objects in the sample). The UV line ratios of three composite spectra of other powerful AGN samples are also shown for comparison: NLRG, radio-loud (RL), and radio-quiet (RQ) QSO1 (Telfer et al. 2002). The NLRG composite (Vernet et al. 2001) was created by combining high S/N Keck spectra of 9 powerful radio galaxies at z ∼ 2.5. The RQ and RL QSO1 composites (Telfer et al. 2002) were created by combining Hubble Space Telescope (HST) spectra of QSO1 at 0.3 ≲ z ≲ 3.6. Most objects (139 out of 184) are at z ≲ 1.5. The RQ and RL composites include spectra of 107 and 77 objects, respectively. The median line ratios of the sample of nitrogen-loud (N-loud) QSO1 from Batra & Baldwin (2014) and the median line ratios of SDSS-III BOSS QSO2 candidates at 2 ≲ z ≲ 4.3 (Alexandroff et al. 2018, Silva et al. 2020) are also shown. The location of the 21 core ERQ in several diagnostic diagrams involving relevant UV lines is shown in Fig. 2.

Fig. 2. Comparison of the core ERQ UV line ratios with N-loud quasars (Batra & Baldwin 2014) and the median ratios of SDSS QSO2 at z >  2 (Group 1: CIV/HeII < 4 and Group 2: CIV/HeII > 4, Silva et al. 2020). HeIIλ1640 is not detected for most core ERQ, therefore lower limits are shown with arrows for the ratios involving this line. The minimum and maximum y- and x-axis values have been selected in all diagrams to span approximately the full range of values shown by most plotted objects. The core ERQ line ratios (measured and lower or upper limits) partially overlap the area of the diagrams covered by N-loud quasars.

3.1.1. Comparison with luminous type 2 AGN

The comparison with NLRG is interesting because these are obscured AGN with quasar-like luminosities (Vernet et al. 2001). The BLR is hidden, and the emission line spectrum originates in the NLR. As shown by Humphrey et al. (2008), the UV and optical emission line ratios of high-z NLRG are best explained by AGN photoionisation of low-density gas (n <  1000 cm⁻³).

The line ratios of the core ERQ are inconsistent with NLRG (Table 1, see also Silva et al. 2020). NLRG show much stronger HeII than other lines and fainter NV. Typically, CIV/HeII < 2, CIII]/HeII < 1, NV/CIV ≲ 1, and NV/HeII ≲ 1 in NLRG (Villar Martín et al. 1997; De Breuck et al. 2000; Humphrey et al. 2008; see also Table 1). Figure 2 shows that in general, core ERQ show CIV/HeII ≳ 3, CIII]/HeII ≳ 1, NV/CIV ≳ 1, and NV/HeII ≳ 3. Star-forming galaxies also show relatively weak HeII, with high CIV/HeII and CIII]/HeII similar to core ERQ. However, stellar photoionisation can be discarded with confidence as a relevant excitation mechanism of the gas on the basis of the huge EW of the emission lines and the strength of NV (e.g. Feltre et al. 2016; Nakajima et al. 2018).

Alexandroff et al. (2013) identified a sample of 145 QSO2 candidates at 2 ≲ z ≲ 4.3 selected from the quasar sample of the SDSS-III BOSS. They have a weak rest-frame UV continuum (typical continuum magnitude of i ∼ 22) and strong CIV and Lyα, with an FWHM ≲ 2000 km s⁻¹. Silva et al. (2020) have recently studied the physical properties and abundances of the ionised gas in these systems based on their location in rest-frame UV diagnostic diagrams and on the comparison with AGN photoionisation model predictions. Compared to NLRG at similar redshifts, the QSO2 are offset to higher NV/HeII, CIV/HeII, and CIII]/HeII.

Silva et al. (2020) have classified the QSO2 into two groups: objects with normal CIV/HeII ratios of < 4 (Group 1), and those with extreme CIV/HeII ratios of > 4 (Group 2). Group 2 QSO2 also has systemically higher NV/HeII ratios for otherwise similar line ratios. To explain such systems as well as the difference with Group 1 QSO2 and high-z NLRG, Silva et al. (2020) propose a combination of high gas density n ≳ 10⁷ cm⁻³ and/or supersolar abundances Z ≳ 4 × Z_⊙ with N/H ≳ 16 × Z_⊙ assuming secondary production of N that results in a quadratic (rather than linear) increase of its N/H abundance. The secondary production of C is also discussed by these authors as a possible explanation for the extreme CIV/HeII ratios seen in some QSO2s. However, in the case of core ERQ, no models can explain the line ratios simultaneously.

Based on the value CIV/HeII > 4, 17 of the 21 core ERQ studied here are confirmed to belong to Group 2 (Fig. 2, tope left panel). The classification of the remaining 4 is not possible based on the CIV/HeII lower limits. However, the high NV/HeII (lower limits in the range ∼4–7) strongly suggest that they also belong to Group 2 (NV/HeII < 3 for Group 1 QSO2, Silva et al. 2020). The same classification appears to apply to core ERQ in general, based on the visual inspection of their spectra (H17).

The four core ERQ studied by Z16 roughly show the lower (J0834+0159 and J2215−0056) and upper values J1232+0912 and J2323−0100) of the range of line ratios spanned by core ERQ (Fig. 2; we recall that only upper limits were measured for HeII). We highlight them with different symbols for guidance and as rough markers of the extreme locations in the diagrams of core ERQ. Most of the remaining core ERQ overlap with them or lie in between.

We next discuss the best-fit AGN photoionisation models within the range explored by Silva et al. (2020). The authors considered n in the range 10^2 − 8 cm⁻³ and the gas metallicity in the range Z = (0.5 − 5.0) × Z_⊙.

J0834+0159 and J2215−0056. These sources have some of the lowest UV line ratios of core ERQ. Based on their location in the diagrams of NV/CIV vs. NV/HeII, NV/CIV vs. CIV/CIII], and CIV/HeII vs. CIV/CIII], the parameters of the best-fit models are listed below (see Fig. 1 in Silva et al. 2020).

• Secondary N and primary C models. Models with n ∼ 10^6 − 8 cm⁻³ and 3–5 Z_⊙ are consistent with the location of the two objects in all diagrams mentioned above. These models are also consistent with NIV]/CIV ≲ 0.3 measured in these two objects.

• No adequate secondary N and secondary C models are found. Only models with n ∼ 10^2 − 4 cm⁻³ and Z ∼ 2 Z_⊙ can reach NV/CIV ∼ 2. However, these fail in other diagrams. For instance, they produce NV/HeII ∼ 1 (compared with the observed NV/HeII > 3). Moreover, they produce CIV/CIII] ≳ 5, compared with the measured CIV/CIII] = 1.7.

Within the range of models explored by Silva et al. (2020), the models favoured for J0834+0159 and J2215−0056 therefore have primary production of C and secondary production of N, n ∼ 10^6 − 8 cm⁻³ and Z ∼ 3 − 5 Z_⊙.

J1232+0912 and J2323−0100. Their line ratios are rather extreme. NV is very strong compared with QSO2: NV/CIV ∼ 2 and NV/HeII ≳ 10, while QSO2 in general show NV/CIV < 1 and NV/HeII < 10. CIV/CIII] ∼ 5 is also at the high end of the range of values spanned by QSO2. None of the best-fit models within the parameter space explored by Silva et al. (2020) can explain all line ratios simultaneously. On the other hand, the high NV ratios and large CIV/HeII strongly suggest high densities n >  10⁶ cm⁻³ and/or high metallicities Z ≳ 4 Z_⊙.

The discussion above suggests that densities n ≳ 10⁶ cm⁻³ and well above solar metallicities are unavoidable in the gas that emits the UV lines in the four core ERQ. Because these ratios span roughly the most extreme values of core ERQ line ratios (Fig. 2), we extend this conclusion to core ERQ in general. We argue below that the UV emission lines in core ERQ have a significant and even dominant contribution from the BLR (see also Alexandroff et al. 2018).

3.1.2. Comparison with QSO1

The relative weakness of HeII is characteristic of QSO1. The UV emission lines in these systems are preferentially emitted in the BLR. They typically show CIV/HeII > 35, CIII]/HeII > 10, and NV/HeII > 15 (Table 1; see also Vanden Berk et al. 2001). Moreover, most core ERQ show higher NV than the QSO1 composites (NV/CIV ∼ 0.45; see also Vanden Berk et al. 2001). The 97 objects in H17 have median NV/CIV = 1.44 and 90% show NV/CIV > 0.5. Sixty-five percent show NV/CIV > 1.0 in the range ∼1.0–2.9.

The enhanced NV emission is reminiscent of N-loud quasars. The strong UV nitrogen lines in these systems have been interpreted in terms of supersolar metallicities in the BLR (e.g. Hamann & Ferland 1993; Dietrich et al. 2003; Baldwin et al. 2003; Nagao et al. 2006; Matsuoka et al. 2011), sometimes as extreme as Z >  10 Z_⊙ (e.g. Baldwin et al. 2003). Batra & Baldwin (2014) analysed the UV spectra of 41 N-loud quasars at z ∼ 2.0–3.5. They reported median NV/CIV = 0.95 and NV/CIV ∼ 1.0–2.4 for half the objects. The inferred metallicities are in the range ∼1–18 Z_⊙ with median 5.5 Z_⊙.

The UV line ratios of the Batra & Baldwin (2014) N-loud QSO1 are shown in the diagnostic diagrams of Fig. 2. The two Telfer et al. (2002) QSO1 composites do not appear because their high ratios relative to HeII place them outside three of the diagrams. Core-ERQ partially overlap the area that is covered by N-loud QSO1 in the diagrams. This relation is likely to apply to the general population of core ERQ because the line ratios suggested by the visual inspection of their spectra (H17) are similar. This suggests that the UV lines in core ERQ are emitted in metal-rich BLR.

Different authors have reported a correlation between quasar L_bol and the BLR metallicity Z, as traced by the nitrogen (N⁺², N⁺³, N⁺⁴) emission line ratios (Hamann & Ferland 1999, Warner et al. 2004, Shemmer et al. 2004, Nagao et al. 2006, Juárez et al. 2009). Because of the high L_bol of core ERQ, high BLR abundances (and the enhancement of the NV emission as a consequence) are similarly expected if they follow the same correlation. According to Xu et al. (2018), Z/Z_⊙ ∼ 7 − 16 Z_⊙ are expected in the BLR of quasars with log(L_bol) ∼ 47.0–48.0 (i.e. similar to core ERQ).

To investigate the origin of the L_bol − Z_BLR relation, Matsuoka et al. (2011) used the optical spectra of 2383 quasars at 2.3 <  z <  3.0 extracted from SDSS. They divided the sample into bins based on M_BH and the Eddington ratio . They then extracted composite spectra for each bin of M_BH and and proceeded to measure the line ratios that are most sensitive to BLR metallicity, which include NV/CIV and NV/HeII. They found that the most lumoninous QSO1 show the highest NV ratios (see Fig. 3). We investigate next how NV/CIV and NV/HeII for core ERQ behave in relation to this QSO1 correlation.

Fig. 3. NV ratios and BLR metallicity. The blue, red, and green symbols correspond to the QSO1 bins in which Matsuoka et al. (2011) distributed their QSO1 sample. Each cross represents the location of one given Lbol bin. The colours indicate ranges of increasing Lbol. The median values of QSO2 (Groups 1 and 2) are also shown for comparison with orange symbols (see text).

Because no reliable M_BH are available for core ERQ (see Sect. 3.3), we instead investigated the behaviour of the line ratios with L_bol. L_bol for the different bins in the Matsuoka et al. (2011) QSO1 sample were obtained using and M_BH in their Table 3. They are in the range 46.34 ≤ log(L_bol) ≤ 47.54. For comparison, 46.98 ≤ log(L_bol) ≤ 47.86 for the core ERQ discussed in this work (Perrotta et al. 2019). We plot the QSO1 bins and the core ERQ ratios in the NV/CIV vs. NV/HeII diagram (Fig. 3). The tight correlation between NV/CIV and NV/HeII for QSO1 at least partially reflects the lack of the Baldwin effect of NVλ1240 in type 1 AGN at z ≲ 5 (Dietrich et al. 2002). While the optical and UV broad and narrow lines of type 1 AGN show an anticorrelation between the REW and the continuum luminosity λL_{1450 Å} (the so-called Baldwin effect, Baldwin 1977), the REW of NVλ1240 remains nearly constant over 6 orders of magnitude in continuum luminosity (Dietrich et al. 2002). The different behaviour of this line has been explained as a consequence of the increasing gas metallicity with AGN luminosity and the secondary production of N; that is, N/O ∝ O/H (Korista et al. 1998).

Core-ERQ show a large scatter, in part due to the non-detection of HeII in most of them. In spite of this, it appears clear that the NV ratios in these systems are among the highest compared to the Matsuoka et al. (2011) QSO1, as we would expect based on their extreme L_bol. This supports the hypothesis that the lines are emitted in metal-rich BLR.

This comparison should be taken with caution. NV ratios are very sensitive to the ionisation level of the gas (Humphrey et al. 2008; Batra & Baldwin 2014). To isolate ionisation effects from those that are due to the gas abundances, it would be necessary to measure other lines such as NIV] and NIII] (Batra & Baldwin 2014), which are in general undetected in the SDSS spectra of the core ERQ. In spite of this limitation, the comparison presented above with QSO1 strongly suggest that the UV lines in core ERQ are emitted by metal-rich BLR (see also H17 and Alexandroff et al. 2018).

3.2. Size and density of the UV emission line region

We argued above that the UV lines are emitted in the BLR of core ERQ. The suppression of the UV/optical continuum, the wingless line profiles, and the intermediate FWHM between QSO1 and QSO2 (see Sect. 3.3) suggest that we observe the outskirts of this region.

The range in density in the BLR comes mainly from the estimated radii and photoionisation theory (e.g. Osterbrock 1989; Ferland et al. 1992). Clouds with densities from ∼10⁹ cm⁻³ (close to n_crit ∼ 3.2 × 10⁹ cm⁻³ of C III]λ1909) to > 10¹³ cm⁻³ are expected. The highest density clouds n >  10¹³ cm⁻³ are continuum radiators (Rees et al. 1989). We therefore assume n ≳ 10⁹ cm⁻³.

The BLR in AGN is photoionised. The main evidence is that the emission-line spectra change in response to changes in the continuum, with lag times corresponding to characteristic radii of the BLR (Peterson 1993). When we therefore assume photoionisation, the radiation pressure likely confines the ionised layer of the illuminated gas (Stern et al. 2014). An implication is that the gas density near the ionisation front varies with distance from the nucleus r as

(1)

where L_i, 45 is the ionising luminosity L_i in units of 10⁴⁵ erg s⁻¹ and r₅₀ is r in units of 50 pc. Knowing n and L_i, 45 we can thus infer distances.

We estimated L_i ∼ 0.35 × L_bol (Stern et al. 2014). If the UV lines we see are preferentially emitted in a region of n ∼ 10⁹ cm⁻³, r_UV ∼ 2.1–6.7 pc or 0.7–2.1 pc for n ∼ 10¹⁰ cm⁻³, depending on the object in the sample of 20 core ERQ in Perrotta et al. (2019) (J1714+4148 is not considered here because L_bol is not available). The radial size of the BLR, r_BLR of z ≲ 0.3 type 1 AGN with log(λL₅₁₀₀)≲45.0 correlates with the rest-frame 5100 Å luminosity λL_λ5100 as (Kaspi et al. 2005). This relation was expanded up to log(λL₅₁₀₀) ∼ 46.0 by Vestergaard (2002). We next constrain r_BLR for the core ERQ, assuming that they follow this correlation.

We inferred λL_λ5100 following log(L_bol) = (4.89 ± 1.66)+(0.91 ± 0.04) × log(λL₅₁₀₀) (Runnoe et al. 2012). If core ERQs follow this correlation, we expect log(λL_λ5100) ∼ 46.1–47.2 for the 20 core ERQ. These are ∼7–45 times higher than the values inferred from the observations (Perrotta et al. 2019). This shows that the optical continuum at 5100 Å is strongly suppressed (H17; Perrotta et al. 2019), as is naturally expected in our proposed scenario. The λL_λ5100 inferred from L_bol imply r_BLR ∼ 0.6–3.1 pc depending on the object, which agrees reasonably well with the r_UV values estimated for n ∼ 10^9 − 10 cm⁻³. Because of their high luminosity, the expected BLR sizes are ∼24–126 times larger than the median size inferred for the low-z type 1 AGN sample (29.2 lt-days, Kaspi et al. 2005). The sizes could be somewhat larger if obscuration makes any of the near-and mid-infrared emission anisotropic because the intrinsic luminosities would be even higher (H17). Core-ERQ are therefore expected to have large BLR because of their extreme luminosities. This may favour the partial visibility at a wider range of inclinations than for less luminous quasars.

3.3. Orientation and partial view of the broad line region

The effect of orientation on the properties of the SED and the emission line spectra of broad line active galaxies (e.g. quasars and broad line radio galaxies, BLRG) has been widely investigated in the literature (e.g. Marziani et al. 2017; Liu et al. 2018). Some authors have proposed that at least a fraction of BLRG are partially obscured misdirected quasars, seen at somewhat larger angles than the quasar population. In this scenario, they are quasars seen through the edge of the obscuring torus (Dennett-Thorpe et al. 2000; Morganti et al. 1997), which therefore does not have a distinct edge (Baker 1997).

The role of orientation in the context of core ERQ has also been discussed (H17). Alexandroff et al. (2018) highlighted its essential role in the model they proposed to explain the UV emission line and continuum spectropolarimetric properties of the core ERQ SDSS J1652+1728 at z = 2.94. The authors describe a polar outflow model seen at an intermediate orientation between type 1 (face-on) and type 2 (edge-on), with the UV lines being produced on spatial scales similar or greater than the scales of the dusty broad emission and broad absorption emission line regions (see also H17). This scenario is in turn inspired by the model proposed by Veilleux et al. (2016) for the nearest quasar MRK231.

We next review and reinforce with new arguments the significant role that orientation plays in determining the observed properties of core ERQ. Based on all arguments combined, we propose that a significant fraction of core ERQ are very luminous but otherwise normal quasars that are seen at an intermediate orientation between type 1 (face-on) and type 2 (edge-on). This results in a partial view of the BLR, with the outskirts being visible, while the inner BLR and the AGN continuum remain hidden. This scenario is consistent with seven aspects that we detail below:

1. The suppression of the UV and optical continuum compared with blue QSO1 (Sect. 3.2, H17)

2. The UV line ratios consistent with metal rich BLR (Sect. 3.1). The continuum source is hidden from view, but the BLR is partially observable.

3. The core ERQ intermediate UV line FWHM values between QSO1 and QSO2 at similar z (Sect. 1) and the lack of broad wings typical of blue QSO1 (H17). This can be explained if the inner BLR is obscured. The inner BLR clouds would be responsible for the broad wings of the UV lines because they are closer to the SMBH and are expected to have higher velocities (gravitational and non-gravitational as a result of, e.g., radiation-driven outflows). In a situation like this, the FWHM of the broad emission lines is no reliable tracer of the BLR kinematics and inadequate for estimating BH masses.

4. The high LOS hydrogen column densities N_H ≳ 10²³ cm⁻² implied by X-ray observations are probably not due to galactic absorption and are naturally explained by dusty obscuring structures of the central AGN (e.g. the classical equatorial optically thick dusty torus and polar dust; Ramos Almeida & Ricci 2017; Goulding et al. 2018). Moreover, the N_H measured for typical type 1 AGN, for which a nearly face-on orientation is expected, are overall at least a factor of ten lower than those measured for core ERQ (e.g. Mateos et al. 2005, Mateos et al. 2010, Corral et al. 2011).

5. The intermediate REW_[OIII] between QSO1 and QSO2. Orientation seems to play a major role (although it is not the only driver) in the variance of the REW_[OIII] in QSO1. REW_[OIII] is used as indicator of the accretion disc inclination (e.g. Risaliti et al. 2011; Bisogni et al. 2017a). Its value increases from face-on to more edge-on systems because the contribution of the AGN continuum and FeII contamination from the BLR decrease.

The [OIII]λ5007 REW of the Perrotta et al. (2019) sample have REW_[OIII] ∼ 42–646 Å with median 177 Å. For comparison, WISE/SDSS selected hyper-luminous (WISSH) QSO1 for which [OIII] could be detected and measured have REW_[OIII]∼ 0.2–80 Å (median 2.8 Å); the dispersion of these values is due at least in part to orientation (e.g. Bischetti et al. 2017; Vietri et al. 2018). Another parameter that influences REW_[OIII] is L_bol. Very luminous systems show lower REW_[OIII] (the Baldwin effect mentioned above). The L_bol of the core ERQ is within the range of the WISSH QSO1 (> 10⁴⁷ erg s⁻¹). If this is underestimated because of obscuration (H17), the large REW for such high intrinsic L_bol further supports a more edge-on inclination than luminous blue QSO1.

6. NV/CIV vs. FWHM_CIV. The general properties (ionisation, density, kinematics, and abundances) of the BLR and the NLR in AGN are highly different, and this is clearly reflected in the BLR and NLR line FWHM and the UV emission line ratios (see comparison between the NLRG and the QSO1 line ratios in Fig. 2). It is natural to expect that from an edge-on to a face-on orientation, the UV emission line properties gradually change from QSO2 type to QSO1 type because the BLR emission becomes more prominent. The UV line FWHM and the NV line ratios are two of the most clearly distinct parameters in type 1 and type 2 AGN spectra. As explained in the previous section, high NV ratios are characteristic of the BLR, but not of the NLR. We therefore expect that FWHM_CIV and NV/CIV both increase from type 2 to type 1 orientation because the BLR gradually becomes more visible and its contribution to the UV spectrum becomes increasingly dominant relative to the NLR.

We assigned the 97 core ERQs in H17 to five groups according to the range of FWHM_CIV: 1000–2000, 2000–3000, 3000–4000, 4000–5000, and 5000–6000 km s⁻¹. We calculated the median value of FWHM_CIV and NV/CIV for all groups and plot them in Fig. 4. The scatter of the whole sample is huge, but when the median values are considered, it is found that NV/CIV increases as the lines become broader. Objects with FWHM <  2000 km s⁻¹ may be similar in orientation to the QSO2 sample of Alexandroff et al. (2013), although they are about ten times more luminous. The large scatter of the data shows that orientation is not the only parameter influencing NV/CIV and FWHM_CIV. Other parameters such as n, metallicity, the ionisation parameter U, and the nature of the kinematics must be involved (as is the case in other AGN). In spite of this complexity, that the median values behave as expected from a gradual change of inclination, adds support to the idea that orientation plays an important role. According to this scenario, groups of increasingly broader CIV have an increasing relative contribution of the BLR to the emission line fluxes, that would result on higher NV/CIV, more typical of QSO1.

Fig. 4. NV/CIV vs. FWHMCIV for the 97 core ERQs in H17. J1714+4148, the four core ERQ from Z16 and the remaining core ERQ in Perrotta et al. (2019) are plotted with the same colour code as in Fig. 1. All other core ERQ from H17 are plotted with small dark red diamonds. The core ERQ are organised into five groups according to the FWHMCIV range of values, indicated with the vertical dashed lines. The median NV/CIV and FWHMCIV in km s−1 of the five groups are shown with large orange circles. The error bars correspond to the standard error on the median values for each group. A clear increase with line width is observed in the NV/CIV ratio.

7. An intermediate orientation may also explain the heavily absorbed Lyα, which suggests that a large amount of absorbing gas is available. The Lyα profile in core ERQ is often absorbed across most or even the entire spectral line profile (Fig. 5; see also H17). It is likely that in at least some objects, the emission at ∼1216 Å flux is not dominated by Lyα but by OV]λλ1214,1218 instead (Humphrey 2019). The Lyα absorbers span a broad range of velocities of thousands of km s⁻¹ that are sometimes both red- and bluesfhited. This is reminiscent of the broad line absorbers (BLA) (also both blue- and/or redshifted depending on the lines) found by Zhou et al. (2019) in a small sample of QSO1. The authors proposed that the highly redshifted BALs arising from neutral hydrogen, helium, and FeII atoms are due to inflowing gas that directly feeds the accretion disc. The authors proposed high inclination angles, intermediate between face-on and edge-on, and a location for the absorbers between the accretion disc and the dusty torus. This gas, possibly mixed with dust, may be the same reservoir that could explain the illuminated phase of a disc-like BLR that is inflated by an accretion disc wind, as proposed by Czerny & Hryniewicz (2011) and Baskin & Laor (2018).

Fig. 5. SDSS spectra of J1031+2903 (top), J1356+0730 (middle), and J1604+5633 (bottom) in the redshifted Lyα+NV region. Lyα is highly absorbed in many core ERQ. The expected Lyα and NVλ1240 central λ are shown based on the CIV and SiIV+OIV] redshifts. Fluxes are plotted in units of ×10−17 erg s−1 cm−2 Å−1. The wavelength is shown in Å.

3.4. [OIII] outflows and the effects of orientation

The unified model proposes that the [OIII]λ5007 emission in AGN originates from regions outside the central obscuring structure. The emission is therefore expected to be independent of orientation. On the other hand, because the AGN continuum is strongly anisotropic, some [OIII] properties do indeed depend on orientation. QSO1 studies have shown that not only REW_[OIII] changes as a function of orientation, but also the prominence and observed kinematics of the outflowing gas in the NLR as traced by the [OIII] profile (e.g. Shen & Ho 2014; Bisogni et al. 2017a; Marziani et al. 2017). The blue broad outflow component decreases both in intensity and in velocity shift relative to the reference wavelength with increasing REW_[OIII]. Shen & Ho (2014), Bisogni et al. (2017a), and Marziani et al. (2017) interpret these results in terms of a gradual change in orientation from face-on (low REW_[OIII]) to more edge-on (high REW_[OIII]) positions.

In an orientation-based scenario, the most extreme kinematics associated with an ionised outflow are expected to be observed in face-on luminous quasars. This means that extreme [OIII] outflows with velocities of several 1000s km s⁻¹ should be detected in QSO1 with L_bol similar to those of core ERQs. In contrast, Perrotta et al. (2019) found that core ERQ exhibit the broadest and more blueshifted [OIII] emission lines ever reported, with outflow velocities W₉₀ about three times higher than those of luminosity-matched blue quasar samples (see also Z16). These extreme and apparently unique [OIII] kinematics contradict the orientation scenario. We revise next whether there is enough evidence to sustain this conclusion.

3.4.1. Fitting [OIII] in luminous QSO1

Perrotta et al. (2019) considered for the comparative QSO1 sample the five WISSH QSOs presented in Bischetti et al. (2017) and the 74 luminous QSO1 in Shen (2016), which have 46.2 ≤ log(L_[OIII]) ≤ 48.2 and 1.5 <  z <  3.5. They re-built the artificial [OIII]λ5007 lines for each QSO1 using the best-fit Gaussian parameters quoted by Shen (2016) and Bischetti et al. (2017). They compared the inferred non-parametric velocities W₉₀ with those of core ERQs. No errors have been estimated. They derived an average W₉₀ = 1550 ± 300 km s⁻¹ for this QSO1 sample.

The first aspect to take into account is that the low REW_[OIII], the complex Hβ+[OIII] blend, and the strong and complex FeII emission in luminous QSO1 makes the determination of [OIII] spectral parameters very challenging. Different works have shown that REW_[OIII] decreases with increasing AGN luminosity, albeit with a very large scatter (the [OIII] Baldwin effect; Brotherton 1996; Dietrich et al. 2002; Shen 2016). That is, if the [OIII] lines are detected at the highest luminosities, it is more difficult (often impossible) to parametrise them because their REW is low and the broad Hβ and the underlying FeII multiplets contaminate them relatively more strongly. We argue that the particular orientation of core ERQ facilitates the detection and parametrisation of the outflows because the contrast of the [OIII] lines is higher with respect to the continuum, the FeII multiplets, and the broad Hβ because these are totally or partially obscured.

We note that first, about 70% of the 18 WISSH QSO1 analysed in Vietri et al. (2018; see also Bischetti et al. 2017) show weak or absent [OIII] emission. This hinders the detection and accurate characterisation of the outflows in most of their targets. Second, about 35% of the QSO1 in Shen (2016) is associated with [OIII] detections with < 3σ significance. Furthermore, the QSO1 in Shen (2016) are characterised by low REW, with an average REW_[OIII] = 14 ± 10 Å (see Fig. 4 in Shen 2016). The outflow velocities must therefore be associated with very large uncertainties.

Recently, Coatman et al. (2019) re-analysed all rest-frame spectra of z >  1.5 QSO1 available from the literature, in order to infer the properties of ionised outflows traced by [OIII] gas with a uniform analysis strategy. The authors provided new independent velocity measurements for the Shen (2016) targets and also reported significant differences in the derived [OIII] fluxes and velocities. They showed that such discrepancies are not due to the quality of the analysed near-infrared spectra, but to the systematic effects associated with the fit. For instance, the templates that were chosen to model the FeII emission (see e.g. Kovačević et al. 2010; Vietri et al. 2018, and references therein), the function(s) adopted to reproduce the BLR emission (e.g. single or multiple Gaussian or Lorentzian functions, or broken power-laws; see e.g. Nagao et al. 2006; Shen et al. 2011), and the shape of the underlying continuum (e.g. Mejiá-Restrepo et al. 2016; Varisco et al. 2018) can significantly affect the analysis results. All these systematic effects are increasingly important for QSO1 with high L_bol, for which we expect more extreme outflows (e.g. King & Pounds 2015; Ishibashi et al. 2018) but also higher emission from the more nuclear regions relative to the [OIII] line flux (e.g. Shen & Ho 2014; Shen 2016).

To further highlight the difficulties in deriving robust [OIII] kinematics of high-luminosity QSO1, we derived new outflow velocity measurements for a subsample of blue QSO1 presented in Shen (2016) using a different method, the multicomponent simultaneous fit analysis (e.g. Brusa et al. 2015). We focused on the 13 out of 19 targets at z >  2 (as the core ERQs) associated with good-quality near-infrared spectra. We refer to Appendix B for details regarding the fit analysis and the main differences with respect to the previous work by Shen (2016) and Coatman et al. (2019). In Table 2 we report the W₈₀ and W₉₀ measurements we derived from the [OIII] profiles, together with those from the analyses of Shen (2016) and Coatman et al. (2019; see also Table B.1.

The three sets of W₈₀ values are also compared in Fig. 6. Our results better agree with those of Coatman et al. (2019) than with those of Shen (2016), who reported line width measurements that are 1.3–3 times narrower. The differences in the analysis results between our work and these two studies highlight the difficulties in deriving robust [OIII] kinematic measurements for very luminous QSO1. This difficulty is also reflected in the large errors of our W₈₀ values. They are larger than those in Coatman et al. (2019). We consider our errors more realistic because they take into account (whenever possible) all systematic effects related to the modelling of continuum, FeII, BLR, and NLR emission (see Appendix B). Non-parametric velocities errors are not available for Shen (2016).

Fig. 6. Comparison of our [OIII]λ5007 W80 with Coatman et al. (2019) and Shen (2016) for 13 luminous QSO1 at z >  2.

The 13 QSO1 discussed above have L_[OIII] ∼ (0.8 − 15.0) × 10⁴⁴ erg s⁻¹. We infer W₉₀ ≥2000 km s⁻¹ for 11 of them (see Table 2). All are above the median W₉₀ = 1568 km s⁻¹ quoted by Perrotta et al. (2019) for the entire QSO1 sample from Shen (2016). Out of the 18 QSO1 within the L_[OIII] range as the subsample fitted by us, they find only three with W₉₀ ≥ 2000 km s⁻¹ (see their Fig. 4). Based on the Shen (2016) fits, they clearly obtain significantly narrower lines than core ERQ for the blue QSO1 comparative sample. With all the caveats mentioned above, it is clear that these kinematic parameters must be affected by very large uncertainties. Unless errors are quoted, any comparison with core ERQ can only be considered tentative.

It is also important to note that the L_[OIII] core ERQ values in Perrotta et al. (2019) are lower limits. They estimated L_[OIII] for each object by multiplying the [OIII] REW by the luminosity at 5007 Å. To obtain the continuum fluxes at this λ, they linearly interpolated the WISE W1 (3.6 μm) and SDSS i or z magnitudes and then converted the derived values at 5007 Å into fluxes. As pointed out by the authors, the luminosities could be on average 6.3 times higher because the median core ERQ SED is suppressed at ∼5000 Å by about 2 mag relative to normal or blue quasars (H17). The correction factor could be even higher (Sect. 3.2).

3.4.2. [OIII] kinematics. Comparison between core ERQ and luminous QSO1

In Fig. 7 we compare the [OIII] kinematics of core ERQs as traced by W₈₀ with other samples of type 1 AGN. In particular, we consider: (1) Perrotta et al. (2019) core ERQ (the four objects from Z16 are highlighted with different symbols as in Fig. 2). W₈₀ were derived rescaling their W₉₀, assuming W₈₀  ∼  0.72 × W₈₀. This is consistent with the value expected for a Gaussian profile and with the W₈₀/W₉₀ ratios in Zakamska & Greene (2014) for type 2 AGN, in Z16 for core ERQs and in Coatman et al. (2019) for QSO1; (2) SDSS J1714+4148 (see Sect. 2); (3) 13 QSO1 from Shen (2016) at z> 2 using the result from our new fits; (4) low z X-ray detected SDSS type 1 AGN at z <  0.8 from Perna et al. (2017a). This sample allows us to highlight the correlation between L_[OIII] and W₈₀ over a range of five dex, as well as its large scatter at fixed luminosity; (5) five WISSH QSO1 from Bischetti et al. (2017). We measured W₈₀ by reconstructing the [OIII] profiles with the Gaussian parameters provided by the authors in their tables; (6) Coatman et al. (2019) sample of luminous QSO1 (45.5 <  log(L_bol) < 49.0) at 1.5 <  z <  4.0. The 13 QSO1 at z >  2 re-analysed in this work were not excluded from the figure (see Fig. 6 and Table 2 for a comparison between their fit results and ours); (7) twenty-two heavily reddened quasar candidates from the UKIDSS-LAS, VISTA VHS, and VIKING imaging surveys at 2.0 ≲ z ≲ 2.6 and 46.0 ≲ log(L_bol) ≲ 48.5 (Temple et al. 2019).

Fig. 7. Location of core ERQ and several type 1 AGN samples in the W80[OIII] vs. L[OIII] plane.

We find that core ERQ fall on the W₈₀ vs. L_[OIII] correlation as defined by different samples of type 1 AGN. This is consistent with the scenario in which more kinematically extreme outflows are triggered by more luminous AGN. Any reasonable correction factor for L_[OIII] due to dust extinction would maintain them within the correlation.

The [OIII] of some very luminous non-obscured QSO1 is as broad as that of core ERQ (Shen 2016; Bischetti et al. 2017; Coatman et al. 2019). Nine out of the 22 heavily reddened quasar candidates (Temple et al. 2019) show [OIII] kinematics similar to some core ERQ (see Fig. 7). The extreme [OIII] kinematics is therefore not exclusive to this object class. When the lower limits on L_[OIII] and the new results on the blue QSO1 sample of Shen (2016) are taken into account, it cannot be claimed that core ERQ show in general more extreme kinematics than blue QSO1 with similar bolometric luminosities.

Extreme [OIII] kinematics is indeed much more frequent in core ERQ (∼100%) than in QSO1, as already pointed out by Perrotta et al. (2019). When all AGN samples are considered, 8 out of the 9 objects with the most extreme kinematics (W₈₀ ≳ 4000 km s⁻¹) are core ERQ. It is unclear whether this is a real intrinsic difference that arises because more extreme outflows are triggered in core ERQ or is rather an artificial effect due to uncertainties that result from all caveats mentioned above and the easier detection and parametrisation of outflows in core ERQ.

Based on the above, we conclude that the extreme [OIII] kinematics of core ERQ do not pose a problem for the orientation-based scenario.

3.4.3. Kinetic power of the [OIII] outflows

High outflow kinetic powers of at least ∼3–5% of L_bol have been inferred for core ERQ (Perrotta et al. 2019). In principle, such powerful winds have the potential to affect the evolution of the host galaxies. For these calculations, sizes R = 1 kpc and densities n = 200 cm⁻³ were assumed. The authors provide uncertainties to account for a possible range n ∼ 100–1000 cm⁻³ and outflow sizes R ∼ 0.5–7 kpc.

The most likely situation is that the ionized outflows are spatially unresolved (Karouzos et al. 2016; Villar Martín et al. 2016; Husemann et al. 2016; Rose et al. 2018; Tadhunter et al. 2018). The integrated spectra used for the outflow measuments are expected to be dominated by the compact NLR where a broad range of densities of up to n ∼ 10⁷ cm⁻³ are likely to exist (see Villar Martín et al. 2016 for a discussion). For a broad distribution of densities, the emission of a certain forbidden line is expected to peak at gas with n similar to its critical density n_crit (e.g. Stern et al. 2014). For the [OIII]λλ4959,5007 lines, n ∼ n_crit = 8.6 × 10⁵ cm⁻³. If as argued in Sect. 3.2 radiation pressure confinement defines the gas density spatial profile, we can infer the size of the [OIII] region for each object using Eq. (1), n_crit and L_ion. The inferred sizes are ∼73–227 pc, depending on the object. With these sizes and n ∼ n_crit, the kinetic powers of the outflows are reduced by a factor of ∼300–977 depending on the object and become ≪1% L_bol. The existence of lower n gas in the outflows is of course possible (e.g. Baron & Netzer 2019; Perna et al. 2017b) and may carry a significant amount of mass. The existing data prevent us from concluding about this, however.

High spatial resolution spectroscopy is essential to resolve the ionized outflows and characterise their sizes, geometry, and density distribution more accurately. In the meantime, the kinetic powers are too uncertain to infer any useful conclusion regarding their potential to affect the evolution of the host galaxies.

4. Summary and conclusions

Core extremely red quasars have been proposed to represent an intermediate evolutionary phase in which a heavily obscured quasar blows out the circumnuclear ISM with very energetic outflows before it becomes an optical quasar. They are also thought to be sites of extreme large-scale (≳1 kpc) outflows that inhibit star formation.

Based on the revision of the general UV and optical emission line properties of core ERQ at z ∼ 2–3, we propose that a high fraction of core ERQ at least are very luminous but otherwise normal quasars that are viewed at an intermediate orientation between type 1 (face-on) and type 2 (edge-on). The general properties are not necessarily indicative of the particular evolutionary phase mentioned above, but are instead strongly determined by orientation and the effects of the high bolometric luminosities.

The calculations and results presented here rest on the assumption that core ERQ follow many correlations that have been identified in other AGN samples. Although this is uncertain, the argumentations suggest that they host large BLR (r_BLR of about several parsec). The highly inclined LOS is such that we observe only part (the outskirts) of the supersolar metallicity BLR. We propose that the inner BLR and the continuum source remain hidden by the torus.

An intermediate orientation is consistent with many properties of core ERQ: the UV line ratios are similar to N-loud QSO1, the large CIVλ1550 and [OIII]λ5007 rest-frame equivalent widths compared with QSO1 of similar L_bol, the intermediate FWHM of the UV lines between QSO1 and QSO2 at similar z, and their wingless profiles. It is also consistent with the large N_H ≳ 10²³ cm⁻² implied by X-ray observations and the heavily absorbed Lyα profiles.

We present a revised comparison of the [OIII] kinematics in 21 core ERQ (20 from Perrotta et al. 2019 and SDSS J1714+4148, based on GTC EMIR near-infrared spectroscopy) with other samples of type 1 AGN. Core-ERQ fall on the W₈₀ vs. L_[OIII] correlation defined by different samples of type 1 AGN. This is consistent with the scenario in which more kinematically extreme outflows are triggered by more luminous AGN.

The extreme [OIII] kinematics of core ERQ are not exclusive of this object class, as demonstrated by our revised spectral fits of luminous non-obscured QSO1 at similar z and other recent works (Bischetti et al. 2017; Coatman et al. 2019; Temple et al. 2019). In contrast to other works (Perrotta et al. 2019), we find that there is no evidence to support that core ERQ in general show more extreme kinematics compared with blue QSO1 of similar bolometric luminosities. The difficulty of detecting and/or the large uncertainties affecting the parametrisation of the [OIII] profiles in the luminous QSO1 comparison sample, and the fact that both L_bol and L_[OIII] of core ERQ are very uncertain (and underestimated), lie at the core of this discrepancy.

Extreme [OIII] kinematics are indeed much more frequent in core ERQ (∼100%), as has been pointed out by other authors (Perrotta et al. 2019). It is unclear whether this is a real intrinsic difference caused by more extreme outflows triggered in core ERQ or is rather an artificial selection effect that is due to the difficulties of detecting and characterising [OIII] outflows in luminous QSO1 in comparison with the much easier detection and parametrisation in core ERQ.

High spatial resolution spectroscopy is essential for resolving the [OIII] outflows and characterising their sizes, geometry, and density distribution more accurately. In the meantime, the kinetic powers are too uncertain to infer any useful conclusion regarding their potential of affecting the evolution of the host galaxies.

As proposed by other authors, core ERQ are ideal systems for detecting and characterising the most extreme AGN outflows. It remains to be decided whether this is a consequence of the blow-out quasar phase or of a particular intermediate orientation.

Acknowledgments

We thank José Acosta Pulido for valuable scientific discussions and the anonymous referee for the thorough revision of the manuscript. Thanks also to Yue Shen for making his QSO1 near-infrared spectra available. Partly based on observations made with the GTC telescope and the EMIR instrument at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias (program GTC14-19A). We thank the GTC staff for their support with the observations. We thank Joel Vernet for providing the high z NLRG composite spectrum. M. V. M. acknowledges support from the Spanish Ministerio de Ciencia, Innovación y Universidades through the grants AYA2015-64346-C2-2-P and PGC2018-094671-B-I00 (MCIU/AEI/FEDER,UE). M. P. is supported by the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. M. P. acknowledges support from the Spanish Ministerio de Economía y Competitividad through the grant ESP2017-83197-P. AH acknowledges FCT Fellowship SFRH/BPD/107919/2015; Support from European Community Programme (FP7/2007-2013) under grant agreement No. PIRSES-GA-2013-612701 (SELGIFS); Support from FCT through national funds (PTDC/FIS-AST/3214/2012 and UID/FIS/04434/2013), and by FEDER through COMPETE (FCOMP-01-0124-FEDER-029170) and COMPETE2020 (POCI-01-0145-FEDER-007672). S. M. acknowledges financial support through grant AYA2016-76730-P (MINECO/FEDER). This research has made use of: (1) the VizieR catalogue access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in Ochsenbein et al. A&AS, 143, 23; (2) data from Sloan Digital Sky Survey. 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. The SDSS Web Site is http://www.sdss.org/; (3) the Cosmology calculator by Wright (2006); (4) the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: UV line ratios of core ERQ

We show in Tables A.1 and A.2 the line ratios for the 16 core ERQ that are not included in Table 1. For all objects, the line fluxes were measured by integrating the flux within the area defined by the line profiles. This method is more accurate than fitting Gaussian profiles when the lines are highly asymmetric due to absorption, for instance. Gaussian fits were also applied when the line of interest was blended with a neighbouring line. This was useful, for instance, for measuring CIII]λ1909 more precisely. This line is often blended with SiIII]λλ1883,1892. In a minority of objects, when the lines are very broad, AlIII]λ1857 can also contaminate this blend. In such cases, different constraints were applied to define a range of possible CIII] fluxes. As an example, when necessary, the FWHM of CIII] was fixed in km s⁻¹ to be the same as CIV and/or SiIV, taking into account the additional broadening because it is a doublet (1907 and 1909 Å). The error bars take into account the range of fits that are acceptable for the different constraints we applied.

Appendix B: Simultaneous multicomponent spectral fit

Obtaining accurate [OIII] kinematic parameters is challenging in luminous QSO1 (see Sect. 3.4.1). To evaluate the typical uncertainties on the derived [OIII] velocities, we reanalysed 14 out of 19 QSO1 in the Shen (2016) sample at z >  2 (as the core ERQs) that are associated with good-quality near-infrared spectra. The remaining objects could not be fitted because the [OIII] lines in the spectra were weak or absent.

We used the simultaneous multicomponent fit technique that is generally adopted in the literature to reduce the degeneracies between FeII, BLR, and NLR emission (e.g. Brusa et al. 2015; Perna et al. 2017a,b). We simultaneously fitted the Hβ+[OIII]+FeII complex together with the Hα+[NII]+[SII] system for the sources at z ∼ 2, and MgII+FeII+FeIII and Hβ+[OIII]+FeII regions for the targets at higher z. The best-fit results are shown in Figs. B.1–B.3 and in Tables 2 and B.1. These figures also show the comparison of the [OIII] best-fit obtained by Shen (2016) with ours (we were unable to provide the same for Coatman et al. 2019 best fit because [OIII] Gaussian parameters are not tabulated in their paper). The spectrum of J2238−0821 is also shown in Fig. B.2. Unlike in Shen (2016), our fits indicate that NLR emission is undetected in this source (consistent with Coatman et al. 2019).

Fig. B.1. Parametrisation of the Hβ+[OIII]+FeII (left and middle panels) and Hα+[NII]+[SII] (right) regions of z ∼ 2 QSO1 from the sample of Shen (2016). The red lines show the best-fit from multiple systemic, outflow, and broad Gaussian profiles to the line features, together with a single power law to the continuum and one (two) template(s) to the FeII features (see text). Each curve with a different colour represents a distinct Gaussian set (or a FeII template) with the same kinematic properties, as labelled in the figure. In the central panel we show for each target the spectrum around the Hβ+[OIII] lines after the subtraction of the continuum. The red curve represents the best-fit retrieval model from MC fits, and the shaded area encompass 68% of the MC sample; the black curve shows the [OIII] best fit obtained by Shen (2016).

B.1. Modelling rest-frame optical spectra

The Hα+[NII]+[SII] and Hβ+[OIII]+FeII systems of J0149+1501, J1421+2241, J1431+0535, J1436+6336, and J1220+0004 are redshifted in the K and H bands, respectively. We simultaneously fitted the continuum and all emission lines to reduce the degeneracies between BLR, NLR, and FeII emission. In short, we fitted a (single) power-law continuum, an optical FeII template (Kovačević et al. 2010), and a combination of Gaussian functions to model BLR and NLR emission lines.

In particular, depending on the complexity of FeII emission, we considered one or two FeII templates (see e.g. Vietri et al. 2018). These were convolved with a Gaussian; the width of this Gaussian, the amplitude normalisations, and velocity offset of the FeII templates were free variables in our fit (see Perna et al. 2017a for further details). Gaussian components were used to reproduce the emission lines. Specifically, we used (1) a systemic narrow (FWHM <  700 km s⁻¹) Gaussian component for the [OIII] doublet and for the Hβ, Hα, [NII], and [SII] doublets associated with unperturbed NLR emission; (2) a broad (FWHM >  2000 km s⁻¹) Gaussian component to fit the BLR Hβ and Hα emission; and (3) an outflow (FWHM >  700 km s⁻¹) Gaussian component for all forbidden and permitted emission lines to reproduce prominent and asymmetric wings associated with outflowing gas. One or more additional sets of Gaussian components were considered to reproduce more complex profiles, for instance, asymmetric BLR emission lines (see e.g. J1421+2241, Fig. B.1) and extremely broad [OIII] profiles (e.g. J0843+0750, Fig. B.2; see also Fig. 11 in Coatman et al. 2019).

Fig. B.2. Parametrisation of the MgII+FeII+FeIII (left) and Hβ+[OIII]+FeII (middle and right panels) of z >  2.5 QSO1 from the sample of Shen (2016). The red lines show the best-fit from multiple systemic, outflow, and broad Gaussian profiles to the line features, together with a two local power laws to the continuum and one (two) template(s) to the FeII features (see text). Each curve with a different colour represent a distinct Gaussian set (or a FeII template) with same kinematic properties, as labeled in the figure. In the right panel we show, for each target, the spectrum around the Hβ+[OIII] lines after the subtraction of the continuum. The red curve represents the best-fit retrieval model from MC fits, and the shaded area encompass 68% of the MC sample; the black curve shows the [OIII] best-fit obtained by Shen (2016). For J2238−0921, we do not detect any [OIII] emission line, consistent with Coatman et al. (2019) and at odds with Shen (2016): the peaks at ∼4900, 5000 and 5150 Å are associated with FeII emission (precisely, the 6S iron line group, according to the ordering in Kovačević et al. 2010).

For each set of Gaussian functions, we constrained the wavelength separation between emission lines within a given set of Gaussian profiles in accordance with atomic physics. This means that we constrained the velocity offset of the outflow Gaussians from the narrow (systemic) components to be the same for all the emission lines. Moreover, each emission line within a given set had the same FWHM. Finally, the relative flux of the two [NII]λλ6548,5583 and [OIII]λλ4959,5007 components was fixed to 2.99, and the [SII] flux ratio was required to be within the range 0.44 <  f(λ6716)/f(λ6731) < 1.42; finally, the Hα/Hβ flux ratio was constrained to be > 2.85 (Osterbrock 1989).

The number of sets used to model the spectra was derived on the basis of the Bayesian information criterion (BIC; Schwarz 1978), which uses differences in χ² that penalise models with more free parameters (see e.g. Harrison et al. 2016; Concas et al. 2019).

B.2. Modelling MgII, FeII, FeIII, Hβ, and [OIII]

For the remaining targets at z >  2.5, the Hα+[NII]+[SII] complex is not covered by available spectra. We therefore fitted Hβ+[OIII]+FeII systems together with the MgII+FeII+FeIII lines. Under the assumption that the Hβ and MgII emission from BLR has the same kinematics (e.g. Shen & Liu 2012; Bisogni et al. 2017b) and that the optical and UV iron emission have the same width and velocity (e.g. Bisogni et al. 2017b), the simultaneous fit technique allows us to reduce the degeneracies between NLR, BLR and iron emission for these spectra as well. We fitted the AGN continuum with two local power laws to the wavelength regions associated with the two systems (i.e. [2600 − 3000] Å and [4500 − 5500] Å). In addition to the optical FeII templates, we also considered one or two FeII+FeIII templates from Popović et al. (2019) to model the iron emission underlying the MgII lines. The UV and optical iron templates were constrained to have the same widths and velocity shifts. We considered a combination of Gaussian sets to best reproduce the BLR and NLR emission (see above) in this case as well.

We note that (1) the MgIIλλ2796,2803 doublet is fitted with a single Gaussian because the lines are unresolved in all spectra: MgII emission is generally dominated by the BLR emission, whose width is much higher than the separation between the two member lines. Moreover, (2) the distinction between systemic and outflow components in our best-fit results must be taken with caution because of the degeneracies in the fit analysis, its dependency on the quality of the spectra (in terms of S/N and resolution), and the unknown intrinsic shape of systemic and outflow emission line profiles (e.g. Liu et al. 2013; Bae & Woo 2016).

B.3. Main differences with respect to Shen (2016) and Coatman et al. (2019)

Emission from the optical continuum and iron lines in Coatman et al. (2019) and Shen (2016) are modelled with a local power-law plus a FeII template from Boroson & Green (1992), using the wavelength regions just outside the Hβ–[OIII] complex. This model is then subtracted before fitting the BLR and NLR emission. In our analysis, all possible emitting contributions are fitted simultaneously, considering the entire wavelength range covered by FeII, Hβ, and [OIII] lines. Moreover, for the sources at z ∼ 2, we used a unique power law to model the continuum in the Hβ–[OIII] and Hα–[NII] regions, in order to be less affected by degeneracies between strong FeII and continuum emission. Finally, we preferred to use the FeII templates by Kovačević et al. (2010), which can allow a more robust separation between [OIII] and FeII emission at ∼5007 Å (see also e.g. Appendix A in Perna et al. 2017a).

Shen (2016) and Coatman et al. (2019) fitted the [OIII] lines with two Gaussian components, one for the core and one for the prominent blueshifted wing. The only tied kinematics are those of the [OIII] and Hβ core components.

B.4. Best-fit analysis and results

The spectral analysis was performed with a python routine using the cap-mpfit procedure (part of the pPXF package by Cappellari 2017) and performing Levenberg-Marquardt least-squares minimisation between the data and the model. The best-fit results are shown in Figs. B.1–B.3 and reported in Tables 2 and B.1.

Fig. B.3. Continuation of Fig. 6.

In order to estimate errors associated with our measurements, we used Monte Carlo (MC) simulations. For each modelled spectrum, we collected the fit results of 100 mock spectra obtained from the best-fit final models (red curves in Figs. B.1–B.3) and adding Gaussian random noise (taking as a reference the noise spectra shown with cyan curves in the figures). The errors reported in Table 2 were derived by taking the range that contains 68.3% of values evaluated from the obtained distributions for each Gaussian parameter or non-parametric measurement. Finally, we note that because line profiles generally are non-Gaussian and much broader than the spectral resolution, we did not correct the observed profiles for instrumental effects and report all values as measured. We also note that our estimated errors take into account (whenever possible) all systematic effects related to the (simultaneous) modelling of continuum, FeII, BLR, and NLR emission. Our best-fit models were obtained for each MC trial by starting from a random initialisation of the model parameters. This allowed us to model the spectrum by taking the possible degeneracies between the different emitting components into account. As a consequence, the errors reported in Tables 2 and B.1 are up to a few orders of magnitude larger than those reported in Coatman et al. (2019). This better reflect the difficulties in such a kinematic analysis.

For the sources for which we could fit both Hα and Hβ, together with a single power-law from ∼4800 to ∼7000 Å for the continuum, we can be quite confident about the contribution of BLR Balmer emission and the underlying continuum; all but J1220+0004 are also associated with weak FeII emission. These sources therefore provide the most robust kinematic measurements for the [OIII] gas. On the other hand, the sources at higher z can be associated with higher systematic effects because of the difficulties in modelling the local continuum underlying the Hβ and MgII lines², and the presence of strong iron emission (at least in the UV part). Overall, the differences in the analysis results between our work, Shen (2016), and Coatman et al. (2019) highlight the difficulties in deriving robust [OIII] kinematic measurements for very luminous QSO1. Below we briefly comment on the best-fit results for each target.

J0149−1501. Hβ and Hα are redshifted at the edges of the H and K bands, respectively, where the filter transmission is very low. As a result, the bluest part of the line profiles cannot be modelled properly. The BLR emission is reproduced with one Gaussian component. Narrow emission can only be well constrained for the [OIII] doublet, which also presents very broad line wings. A prominent outflow component is required to reproduce the profiles of all emission lines (although the strong degeneracies do not allow us to constrain this component well for the Balmer lines). Our reconstructed [OIII] profiles present more pronounced blue and red wings than those in Shen (2016) and Coatman et al. (2019), although all non-parametric velocities are consistent within 1σ.

J1421+2241. The BLR is modelled with two Gaussian profiles; the NLR is fitted with a systemic and an outflow component. The blue and red wings of the [OIII] lines are both blended with FeII and Hβ emission. The [OIII] reconstructed profile is very similar to the profiles in Shen (2016) and Coatman et al. (2019).

J1431+0535. The Hβ line is severely affected by sky-line residuals; the broad Hα from the BLR is strongly asymmetric. We modelled the BLR emission with two Gaussian components. Narrow emission can only be constrained well for the [OIII] doublet, which also presents very broad line wings. Our reconstructed OIII profiles present more pronounced blue and red wings than those in Shen (2016) and Coatman et al. (2019), although the derived non-parametric velocities are consistent within 1σ with theirs.

J1436+6336. The BLR lines are modelled with two Gaussian profiles; the NLR is fitted with a systemic and an outflow component. The blue and red wings of the [OIII] lines are both mildly blended with FeII and Hβ. Our reconstructed [OIII] profiles present more pronounced blue and red wings than those in Shen (2016) and Coatman et al. (2019), although all non-parametric velocities are consistent within 2σ.

J1220+0004. The Hα+[NII] complex is severely affected by telluric features; a significant portion of the profile is therefore masked during the spectral fitting. The BLR line profile is reconstructed on the basis of the observed Hβ profile, with two BLR Gaussian sets and one Lorentzian profile. The NLR emission is fitted with a systemic component plus two very broad outflow Gaussians. Two FeII templates are also required to reproduce the Hβ+[OIII]+FeII complex. Our reconstructed [OIII] profiles present more pronounced blue and red wings than those in Shen (2016) and Coatman et al. (2019). Under the assumption that Hα and Hβ have very similar profiles, the flux excess around ∼5000 Å cannot be associated with a Hβ red wing because of the absence of a similar feature in Hα. On the other hand, FeII templates are not capable to reproduce the broad emission close to the [OIII] doublet. These arguments allow us to prefer a best-fit with very broad [OIII] lines. We note, however, that strong degeneracies might still be present in our fit, and that systematic errors might be significant for this target.

J2238−0921. This target shows very complex profiles in the Hβ+[OIII]+FeII and MgII+FeII+FeIII regions. We fit the MgII and Hβ BLR emission with two Gaussian components. The iron emission is modelled with two templates. In the optical region, to reproduce the strong peak at ∼5250 Å, we include in the fit the FeIII and FeIV emission lines at 5270 and 5236 Å (but a strong degeneracy between these lines and FeII features is present at these wavelengths). No NLR emission associated with [OIII] is detected for this source, consistent with Coatman et al. (2019, see the detailed discussion in their paper) and at odds with Shen (2016). This source can be considered as a weak [OIII] QSO1, with a CIV line blueshift of ∼5500 km s⁻¹ and REW_CIV ∼ 10 Å (see Vietri et al. 2018).

J0250−0757. The BLR emission can be well reproduced with a single Gaussian profile; [OIII] is modelled with two Gaussians; kinematic results are consistent with those of Shen (2016) and Coatman et al. (2019).

J0844−0503. The [OIII] profiles have similar characteristics to the core ERQs. We modelled the Hβ and MgII BLR emission with a single Gaussian; the NRL emission is fitted with two very broad outflow Gaussians. Our reconstructed [OIII] profiles present more pronounced blue and red wings than those in Shen (2016) and Coatman et al. (2019).

J0942−0422. The BLR emission can be well reproduced with a single Gaussian profile; the [OIII] is modelled with three Gaussians. Our reconstructed [OIII] profiles present more pronounced red wings than those in Shen (2016) and Coatman et al. (2019).

J0304−0008. The BLR emission can be well reproduced with a single Gaussian profile; the [OIII] is modelled with three Gaussians. Our reconstructed [OIII] profiles present more pronounced wings than those in Shen (2016), consistent with Coatman et al. (2019).

J0843+0750. The BLR emission can be well reproduced with a single Gaussian profile; the [OIII] is modelled with two Gaussians. Our reconstructed [OIII] profiles present more pronounced red wings than those in Shen (2016), consistent with Coatman et al. (2019).

J1019+0254. The BLR emission can be well reproduced with a single Gaussian profile; the [OIII] is modelled with two Gaussians. Our reconstructed [OIII] profiles present more pronounced wings than those in Shen (2016), consistent with Coatman et al. (2019)

J0259+0011. The BLR emission can be well reproduced with a single Gaussian profile; the [OIII] is modelled with two Gaussians. Our reconstructed [OIII] profiles present more pronounced wings than those in Shen (2016), consistent with Coatman et al. (2019)

J1034+0358. The BLR emission can be well reproduced with a single Gaussian profile; the [OIII] is modelled with two Gaussians. Our reconstructed [OIII] profiles present more pronounced wings than those in Shen (2016) and Coatman et al. (2019).

All Tables

All Figures

Fig. 1. GTC EMIR spectrum of J1714+4148 in the Hβ and [OIII] spectral region (black lines). The flux is in units of ×10−18 erg s−1 cm−2 Å−1. Different fits (red lines) were attempted by applying a variety of kinematic constraints to [OIII] and/or Hβ (fits 2–7) or by leaving the kinematic parameters free (fit 1). For each fit, the constraints on the observed FWHM in Å are shown. It is also specified whether Hβ was forced to have the same redshift as [OIII]λ5007 (λ forced, fits 2 and 2b). The reduced chi-squared χ2 are also shown, and the best-fit profiles are presented underneath each panel.

In the text

Fig. 2. Comparison of the core ERQ UV line ratios with N-loud quasars (Batra & Baldwin 2014) and the median ratios of SDSS QSO2 at z >  2 (Group 1: CIV/HeII < 4 and Group 2: CIV/HeII > 4, Silva et al. 2020). HeIIλ1640 is not detected for most core ERQ, therefore lower limits are shown with arrows for the ratios involving this line. The minimum and maximum y- and x-axis values have been selected in all diagrams to span approximately the full range of values shown by most plotted objects. The core ERQ line ratios (measured and lower or upper limits) partially overlap the area of the diagrams covered by N-loud quasars.

In the text

	Fig. 3. NV ratios and BLR metallicity. The blue, red, and green symbols correspond to the QSO1 bins in which Matsuoka et al. (2011) distributed their QSO1 sample. Each cross represents the location of one given Lbol bin. The colours indicate ranges of increasing Lbol. The median values of QSO2 (Groups 1 and 2) are also shown for comparison with orange symbols (see text).
In the text

Fig. 4. NV/CIV vs. FWHMCIV for the 97 core ERQs in H17. J1714+4148, the four core ERQ from Z16 and the remaining core ERQ in Perrotta et al. (2019) are plotted with the same colour code as in Fig. 1. All other core ERQ from H17 are plotted with small dark red diamonds. The core ERQ are organised into five groups according to the FWHMCIV range of values, indicated with the vertical dashed lines. The median NV/CIV and FWHMCIV in km s−1 of the five groups are shown with large orange circles. The error bars correspond to the standard error on the median values for each group. A clear increase with line width is observed in the NV/CIV ratio.

In the text

	Fig. 5. SDSS spectra of J1031+2903 (top), J1356+0730 (middle), and J1604+5633 (bottom) in the redshifted Lyα+NV region. Lyα is highly absorbed in many core ERQ. The expected Lyα and NVλ1240 central λ are shown based on the CIV and SiIV+OIV] redshifts. Fluxes are plotted in units of ×10−17 erg s−1 cm−2 Å−1. The wavelength is shown in Å.
In the text

	Fig. 6. Comparison of our [OIII]λ5007 W80 with Coatman et al. (2019) and Shen (2016) for 13 luminous QSO1 at z >  2.
In the text

	Fig. 7. Location of core ERQ and several type 1 AGN samples in the W80[OIII] vs. L[OIII] plane.
In the text

Fig. B.1. Parametrisation of the Hβ+[OIII]+FeII (left and middle panels) and Hα+[NII]+[SII] (right) regions of z ∼ 2 QSO1 from the sample of Shen (2016). The red lines show the best-fit from multiple systemic, outflow, and broad Gaussian profiles to the line features, together with a single power law to the continuum and one (two) template(s) to the FeII features (see text). Each curve with a different colour represents a distinct Gaussian set (or a FeII template) with the same kinematic properties, as labelled in the figure. In the central panel we show for each target the spectrum around the Hβ+[OIII] lines after the subtraction of the continuum. The red curve represents the best-fit retrieval model from MC fits, and the shaded area encompass 68% of the MC sample; the black curve shows the [OIII] best fit obtained by Shen (2016).

In the text

Fig. B.2. Parametrisation of the MgII+FeII+FeIII (left) and Hβ+[OIII]+FeII (middle and right panels) of z >  2.5 QSO1 from the sample of Shen (2016). The red lines show the best-fit from multiple systemic, outflow, and broad Gaussian profiles to the line features, together with a two local power laws to the continuum and one (two) template(s) to the FeII features (see text). Each curve with a different colour represent a distinct Gaussian set (or a FeII template) with same kinematic properties, as labeled in the figure. In the right panel we show, for each target, the spectrum around the Hβ+[OIII] lines after the subtraction of the continuum. The red curve represents the best-fit retrieval model from MC fits, and the shaded area encompass 68% of the MC sample; the black curve shows the [OIII] best-fit obtained by Shen (2016). For J2238−0921, we do not detect any [OIII] emission line, consistent with Coatman et al. (2019) and at odds with Shen (2016): the peaks at ∼4900, 5000 and 5150 Å are associated with FeII emission (precisely, the 6S iron line group, according to the ordering in Kovačević et al. 2010).

In the text

	Fig. B.3. Continuation of Fig. 6.
In the text