Exploring the Diversity of Type 1 Active Galactic Nuclei Identified in SDSS-IV/SPIDERS

We present a statistical analysis of the optical properties of an X-ray selected Type 1 AGN sample, using high signal-to-noise ratio (S/N>20) spectra of the counterparts of the ROSAT/2RXS sources in the footprint of the SDSS-IV/SPIDERS (Spectroscopic IDentification of eROSITA Sources) programme. The final sample contains 2100 sources. It significantly extends the redshift and luminosity ranges (z ∼ 0.01−0.80 and L0.1−2.4 keV ∼ 2.0×1041−1.0×1046 erg s−1) used so far in this kind of analysis. By means of a Principal Component Analysis, we derive Eigenvector (EV) 1 and 2 in an eleven dimensional optical and X-ray parameter space, which are consistent with previous results. The validity of the correlations of the Eddington ratio L/LEdd with EV1 and the black hole mass with EV2 are strongly confirmed These results imply that L/LEdd and black hole mass are related to the diversity of the optical properties of Type 1 AGN. Investigating the relation of the width and asymmetry of Hβ and the relative strength of the iron emission rFeII, we show that our analysis supports the presence of a distinct kinematic region: the Very Broad Line Region. Furthermore, comparing sources with a red-asymmetric broad Hβ emission line to sources for which it is blue-asymmetric, we find an intriguing difference in the correlation of the FeII and the continuum emission strengths. We show that this contrasting behaviour is consistent with a flattened, stratified model of the Broad Line Region, in which the FeII emitting region is shielded from the central source.


INTRODUCTION
The exploration of the multi-wavelength emission of Active Galactic Nuclei (AGN), including their spectral features, has motivated the development of unification schemes for this class of persistent and extremely luminous objects (e.g., Antonucci 1993;Urry & Padovani 1995;Netzer 2015;Padovani et al. 2017, see also Elitzur 2012 for some caveats about such unification schemes). Simplifying the complex and historically developed AGN zoo is an appealing idea, since most of the identified species share common spectral properties E-mail: jwolf@mpe.mpg.de over large ranges of luminosity and redshift. In this context, two main sub-classes have been defined based on their optical emission line properties: Type 1, which have both broad and narrow optical/UV emission lines, and Type 2, which show only narrow lines. These two classes are unified via the hypothesis that the broad lines in Type 2 sources are obscured (e.g., Netzer 2008;Hickox & Alexander 2018). There is nonetheless substantial variety in the properties within the Type 1 population. The broad lines encode information about the geometry and the kinematics of the inner regions of the system : the Broad Line Region (BLR, e.g., Rees et al. 1989;Peterson 2006;Gaskell 2009). Photoionised by the continuum emission of the central source, the BLR emits transition lines, which are expected to be broadened by the motion of the gas (e.g., Baldwin et al. 1995;King 2016). The virial broadening of the emission lines can, under certain assumptions, be used to estimate the mass of central supermassive black holes (SMBH). This mass estimation method is calibrated using estimates of the BLR radius from reverberation mapping. (Bahcall et al. 1972;Blandford & McKee 1982;Peterson 1993;Gebhardt et al. 2000;Kaspi et al. 2005;Bentz & Katz 2015;Shen et al. 2015). The response delay of the BLR emission to variations in the ionising continuum emission from the central source (i.e., the accretion disk) is used to determine the distance between the BLR and the nucleus, assuming a Keplerian velocity field.
Broad, low-ionisation lines, such as Hβ (rest-frame wavelength at 4862Å) and the MgIIλ2800Å doublet 1 , are generally thought to yield more reliable virial broadening estimators than, for instance, the high ionisation transition CIVλ1549Å (e.g., Trakhtenbrot & Netzer 2012;Mejia-Restrepo et al. 2016;Mejía-Restrepo et al. 2018). This is in part due to the presence of prominent blueshifted components in the CIVλ1549Å line profiles of strong accretors, which may not be related to the orbital motion of the gas and would affect the inferred velocities (e.g., Richards et al. 2002;Zamanov et al. 2002;Marziani et al. 2017;Sulentic et al. 2017).
However, Hβ lines possess profiles which are difficult to reconcile with simple Keplerian motion around the central black hole, and could significantly affect the black hole mass estimates (Negrete et al. 2018). Accurate black hole mass measurements are important for the development of galaxy evolution models. It is therefore critical to understand the mechanisms responsible for the observed broad line shapes.
The debate on the actual BLR kinematics, physical processes and geometry has flourished through the past decades; the variety of the often asymmetric (sometimes doublepeaked) broad lines and their variability has lead to an abundant diversity of competing models. It was not until recently that the kinematics of this region could be resolved for the first time, favouring a flattened BLR as an orbiting extension of the accretion disk (Sturm et al. 2018).
A flattened BLR geometry was initially supported by the relation between the width of Hβ and the radio coreto-lobe ratio. (Wills & Browne 1986;Jarvis & McLure 2006;Brotherton et al. 2015) and by the double shoulders and peaks observed in certain broad emission lines (Oke 1987;Eracleous & Halpern 1994Storchi-Bergmann et al. 2016). The double peaked profile variability has been shown not to correlate with continuum changes. The peak shifts are therefore not a reverberation effect (Wanders & Peterson 1996). The strong variability of the double-peaked emission profiles has also been discussed in the light of different disk models by Lewis et al. (2010), which can account for some of the observed line profiles.
A potentially useful approach to constrain and explore the inner structure of AGN is to acquire an understanding of their spectral diversity. For instance, in their seminal paper Boroson & Green (1992) performed a Principal Component Analysis (PCA) of a set of optical, radio and X-ray features of 87 Type 1 AGN. Through the orientation of the first principal component (Eigenvector 1,EV1) in their parameter space, they established an anti-correlation of the equivalent width of the narrow [OIII] lines at 4959Å and 5007Å with the relative strength of the iron emission r FeII = F(FeIIλ4570)/F(Hβ), and the full width at half maximum (FWHM) of the Hβ line as central markers of diversity in their sample. A large range of subsequent studies have investigated the EV1 correlation planes (e.g., Sulentic et al. 2000a,b;Marziani et al. 2001;Shang et al. 2003;Grupe 2004;Kuraszkiewicz et al. 2009;Mao et al. 2009;Tang et al. 2012). In particular, Sulentic et al. (2000a,b) and Marziani et al. (2001) established the foundations of the four-dimensional EV1 (4DE1) formalism, including FWHM Hβ and r FeII as two of the main correlates of EV1 in Boroson & Green (1992). These two quantities are respectively related to the black hole mass (since FWHM Hβ is used as virial broadening estimator) and the Eddington ratio (e.g., Marziani et al. 2003b;Netzer & Trakhtenbrot 2007;Sun & Shen 2015).
The 4DE1 parameter space was extended with the soft X-ray photon index Γ soft (Boller et al. 1996) and the centroid shifts from rest-frame wavelength of CIVλ1549 at 50% fractional intensity (Sulentic et al. 2007b). The careful exploration of this parameter-space led to a proposed two population paradigm in the low-redshift universe, with an empirical separation at FWHM Hβ ≈ 4000 km s −1 (Sulentic et al. 2000a,b, see section 5 of Marziani et al. 2018 and references therein for a full review of the evidence supporting the population A/B classification.). Zamfir et al. (2010, henceforth Z10) studied the characteristics of the Hβ emission of 477 optically selected AGN (z < 0.7) in the optical plane of 4DE1 (i.e., FWHM Hβ vs. r FeII ). They reported that the line shapes distinctively changed along the sequence, the asymmetry and shifts scaling with both of the optical dimensions.
In this work, we explore the EV1 correlation space exploiting high signal-to-noise ratio optical spectra for the largest sample of X-ray selected AGN used to date for this type of analysis.
The X-ray selection allows to construct a less biased sample of AGN, since the host X-ray related processes are believed to be particularly weak compared to high-energetic AGN emission. At the same time X-rays are able to penetrate the host galaxies and in general large column densities, thus making the sample more complete than with other AGN selection techniques (e.g., Singh 2013; Padovani et al. 2017). The AGN in our sample were detected by ROSAT (Boller et al. 2016;Voges et al. 1999), while their multi-wavelength counterparts were determined in Salvato et al. (2018) and spectroscopically followed-up in the SDSS-IV/SPIDERS program (Dwelly et al. 2017). These sources are exclusively Type 1 AGN and were presented in Coffey et al. (2019, henceforth C19) where sources were selected according to the width of their broad lines.
The final sample contains 2100 sources, which is a factor of 4-18 larger than samples used in previous studies of this type. It significantly extends the probed redshift and luminosity ranges (z ∼ 0.01 − 0.80 and L 0.1−2.4 keV ∼ 2.0 × 10 41 − 1.0 × 10 46 erg s −1 ). We perform a PCA on the optical and X-ray features of these sources, which allows us to determine the markers of spectral diversity in our sample. We specifically investigate the role played by the Eddington ratio and the black hole mass for the total variance in our data. The asymmetry of the broad Hβ contributes strongly to the orientation of EV1 through the chosen spectral parameter space. We thus follow in the footsteps of Z10 and connect our global statistical results to the shape of the broad Hβ emission line, which is discussed in the context of BLR stratification and low-ionisation outflows.
The paper is organised as follows: In Sec. 2 we briefly present the SPIDERS AGN value added catalogue of C19 and how we created our subs-ample. We also outline how the additional properties, which were not provided by the catalogue, were determined. In Sec. 3 we describe the distribution of the emission line asymmetry index. The core results of the direct correlation analysis and PCA on our sample are presented in Sec. 4. We demonstrate how the derived physical parameters Eddington ratio and central black hole mass scale with the obtained principal components in Sec. 5. The distribution of our sample in the EV1 plane and the possible presence of a distinct kinematic BLR region is investigated in Sec. 6. We divided the sample in two distinct subsets according to the sign of their Hβ asymmetry indices. From the comparison of parameter correlations in the two sub-samples we detect an interesting contrast in the scaling relation between source luminosity and the equivalent width of the iron emission. This is discussed in Sec. 7. Conclusions are presented in Sec. 8 Throughout this work, we adopt a standard cosmology: Ω M = 0.3, Ω Λ = 0.7 and H 0 = 70 km s −1 Mpc −1 .

DATA
The sample of Type 1 AGN analysed in this work is extracted from the SDSS-IV/SPIDERS AGN catalogue, presented in C19. The original catalogue compiles spectral information through SDSS DR14 2 (Abolfathi et al. 2018), for 7344 Type 1 AGN detected in the Second ROSAT All-Sky Survey ROSAT/2RXS, Boller et al. 2016) and 1157 Type 1 AGN in the XMM-Newton Slew 1 catalogue (XMMSL1, Saxton et al. 2008).
The multi-wavelength counterparts of these X-ray sources were determined via the Bayesian-based code NWAY by (Salvato et al. 2018), using a prior (Dwelly et al. 2017) based on AllWISE data (Cutri & et al. 2013).
Optical spectra of the counterparts were obtained from two different spectrographs, SDSS and BOSS. These instruments cover different optical wavelength ranges: 3800 − 9200Å for SDSS (counting both channels, with a spectral resolution ranging from 1850 to 2200) and 3600 − 10400Å for BOSS (spectral resolution of 1560-2270 in the blue channel, 1850-2650 in the red channel). A technical summary of the SDSS-IV survey is provided by Blanton et al. (2017). The instrument is presented in Gunn et al. (2006). Finally, a detailed description of the SDSS and BOSS spectrographs is given by Smee et al. (2013).
For the sources in the C19 catalogue, the authors fit the spectroscopy following a procedure detailed described in C19. Here we report only the key features.
In C19, the spectral regions centered around the H β and MgII lines (4420−5500Å and 2450−3050Å, respectively) were 2 SDSS DR14 VACs fit with a multi-component continuum model and a series of Gaussian functions. The present work will specifically focus on the broad Hβ emission line and the narrow [OIII]λ4959Å and [OIII]λ5007Å forbidden transition lines. C19's fitting algorithm used up to four Gaussian functions to fit the Hβ line. The best-fit model (and thus the required number of Gaussian components) was automatically selected on the basis of the Bayesian Information Criterion (BIC). Of the four Gaussians, one accounts for the narrow core, while the remaining potential three components of Hβ are defined as broad if they fulfill the criterion: FWHM Hβ > 800 km s −1 . This relatively complex model allows one to trace not only typical broad bimodal profiles above ∼ 1000km s −1 , but also distinct broader components (see C19, Section 8.1, and see Section 6 of this work). For instance Marziani et al. (2010) propose a broad line decomposition into broad, very broad and blue components, which the model defined in C19 may individually trace.
Up to two Gaussians were used for each of the [OIII] lines: one fitting the narrow core and one tracing shifted wings. The continuum model was constructed from a powerlaw, a host galaxy component and an FeII template. The iron emission and the galaxy contribution are obtained respectively from a normalized Narrow-Line Seyfert 1 I Zw 1 template (Boroson & Green 1992) and an early type galaxy template 3 . Morphological studies (e.g., Grogin et al. 2005) have shown that AGN are mainly found in bulge dominated galaxies. Bulges are predominantly home to old star populations. C19 thus justify the use of an early-type galaxy template for the spectral host contribution by the fact that SDSS fibers sample only the central regions of targeted galaxies. The early-type approximation for the host contribution was also used in previous work (e.g. Calderone et al. 2017). The template has Hβ absorption features. Comparing a subsample with strong host contribution to a sub-sample with low host contribution at z<0.2, we found no significant impact of these features on Hβ line shape diagnostics.
The fitting parameters are listed in the catalogue of C19 along with monochromatic continuum and X-ray luminosities and derived parameters such as bolometric luminosities, estimates for the black hole masses and Eddington ratios.
The single-epoch mass estimation method was used for the derivation of the black hole masses (Vestergaard 2002;McLure & Jarvis 2002;Vestergaard & Peterson 2006;Assef et al. 2011;Liu 2012 andShen 2013 for a review). The masses stem from the Assef et al. (2011) calibration, which builds on the FWHM of Hβ and the BLR radiusluminosity relation from Bentz et al. (2009). For this paper, we have reconstructed the best-fit models from the parameters listed in C19. Two examples of Hβ fit reconstructions are shown in Appendix A of the present work. Details of the broad Hβ line decomposition are also presented.
In the following section we describe the sample construction for our statistical analysis.

Sample construction and selection pipeline
We are primarily interested in Type 1 AGN spectral properties which can be associated with the physics of the BLR Figure 1. Distribution of the median S/N ratios per resolution element of the 2100 SDSS DR14 sources in our final sample(white). Also presented are sub-samples of the full SPIDERS sample in C19 after two consecutive quality cuts: the presence of a 0.1 − 2.4 keV flux and the redshift restriction z < 0.8. and the Narrow Line Region (NLR). We thus assembled a parameter subset inspired by the classical quasar PCA papers by Boroson & Green (1992) and Grupe (2004). They included emission properties of Hβ, FeII, [OIII] line at 5007Å 4 and HeIIλ1640Å, as well as X-ray and optical monochromatic luminosities to trace the emission of the hot electron corona and the accretion disk.
The inclusion of X-ray luminosities forces a restriction to just one of the two X-ray surveys compiled in the C19 catalogue, since 2RXS and XMMSL1 soft X-ray fluxes are measured in different bands, at different times, from two different satellites/instruments and different assumptions used for computing the flux (see Saxton et al. 2008;Boller et al. 2016). We confined the analysis to sources detected in 2RXS only (7344 sources), as they were six times more numerous than XMMSL1 sources in the SPIDERS AGN sample. More specifically we retained 2RXS sources, for which the X-ray flux in the 0.1 − 2.4 keV band is provided in C19 (7154 sources). The missing 190 sources have no X-ray flux measurement listed. They correspond to sources with low photon count-rates.
Since we use the monochromatic luminosity at 5100Å to trace the continuum emission, we limit the sample to z=0.8 so that the rest-frame 5100Å are covered by the spectral window.
The redshift cut ensures the presence of fits of the Hβ centered region. These constraints yield 5926 sources.
We further limit our analysis to sources with good SDSS signal-to-noise ratio per resolution element in order to im- prove the significance of our results. Our analysis relies on the accurate measurement of line shapes, and low-quality spectra might strongly affect the fits. For this reason we consider only sources with median S/N ≥ 20. This critical restriction drops our total source count to 2124 objects. The SDSS signal-to-noise ratio distribution of the original sample and its consecutive cuts are shown in Fig. 1.
Finally, 24 sources were excluded due to a lack of black hole mass estimates, extreme uncertainties on the fit parameters of Hβ or absent [OIII]. A visual inspection of these sources' spectra revealed that these issues might be linked to strong continuum emission.
For the 2100 sources, the following subset of spectral properties was compiled (unless specified otherwise, all wavelengths, energies, equivalent widths are defined in the rest-frame):  It is worth to highlight the following points: • Sources for which the parameters relative to the FeII emission were not listed in C19, typically correspond to AGN with very weak FeII emission. We manually set their r FeII parameters and iron equivalent widths to zero. We thus implicitly assume that their iron emission is too weak to be disentangled from the AGN continuum emission.
• All the spectral properties of our defined sub-set are available in C19 except for the equivalent widths and asymmetry indices which were derived as presented in Sec. 2.2 and Appendix A.
• The L 5100Å provided in C19 includes the host contribution. However, in the present analysis we aim at tracing the accretion disk emission with L 5100Å . We thus derived the monochromatic luminosity from the reconstructed power-law model at rest-frame, removing the FeII and host contribution. The uncertainties are obtained from the errors on the normalization of the slope. The considerable degeneracy between the host, iron and power-law components at lower redshift (z < 0.2) due to the stronger host contribution makes it challenging to correctly perform an AGNhost decomposition. It is the principal limitation of deriving monochromatic luminosities directly from the fitted powerlaw model. The relative contributions of the FeII complex and the host galaxy to the monochromatic flux a 5100Å are illustrated in Fig. 2. While the FeII contribution at 5100Å to the total monochromatic flux is marginal over the complete redshift range, the host galaxy strongly contributes to the monochromatic flux at lower redshifts.
• The X-ray fluxes were obtained by C19, following Dwelly et al. (2017), from an absorbed power-law (Γ = 2.4). In addition to these classical fluxes, C19 provide Eddington bias corrected 2RXS fluxes using a Bayesian prior (Kraft et al. 1991;Laird et al. 2009;Georgakakis & Nandra 2011). The Bayesian flux estimates reach unrealistic low values for sources with low count rates. C19 allow up to a factor ten difference with the classical fluxes. Sources for which the Bayesian fluxes do not meet this requirement only have classical flux measurements. We chose to use the uncorrected X-ray flux measurements in this work, since using Bayesian fluxes would reduce the size of our final sample by a factor of two. We note that a test run of the statistical pipeline presented in this work on a smaller sample with Bayesian soft X-ray luminosities yielded similar results to those presented here for the entire sample using the nominal fluxes. Fig. 3 presents the soft X-ray luminosity-redshift distribution of our sources, and compares our sample to that of Grupe (2004). Our sample spans a range of soft X-ray luminosities from 1.9 × 10 41 erg s −1 to 9.9 × 10 45 erg s −1 , with redshifts up to 0.8. Compared to Grupe (2004), we sample to lower luminosities and higher redshifts. The larger sample will in particular allow placement of more stringent constrains on the relation of EV1 and EV2 to physical parameters (see Fig. 11).

Measuring asymmetry in emission lines
This section presents the derivation of a key emission line diagnostic: the asymmetry index.
Emission lines in AGN show asymmetries which encode precious information on the geometry and kinematics of the emitting region. In the case of broad emission lines, asymme- Figure 3. The observed soft X-ray luminosity-redshift for our 2100 2RXS sources. The red triangles are the 110 AGN used in Grupe (2004), which all lie at the brighter end of our sample for equivalent redshifts. Our sample significantly extends both the redshift and the luminosity ranges.
Sample cuts (Initial SPIDERS AGN sample: 7790 sources) remaining sources 2RXS flux 7154 z < 0.8 5926 S/N ≥ 20 2124 Robust asymmetry index/black hole masses 2100 Table 1. Summary of the sample construction pipeline try arises from the relative displacement of the profile's peak and base component's centroid wavelengths. The peak wavelength of the profile is in practice measured at a high fraction of the broad component's intensity (without the narrow core) and is related to the classical broad component (Brotherton et al. 1994;Popovic et al. 2002;Adhikari et al. 2016). The emission line's profile base is, as its name indicates, measured at lower fractional intensities and is expected to trace a distinct, red-or blue shifted, sometimes broader component. In Marziani et al. (2010) and subsequent works of this group, the authors distinguish blueshifted and broad redshifted components of the full profile ('BLUE' and Very Broad Component 'VBC'; for a review see Section 6.2 in Marziani et al. 2018). The shifted VBC is believed to originate from the presence of a distinct emitting region in the BLR: the Very Broad Line Region (VBLR). The VBLR is expected to be composed of highly ionised gas situated even closer to the central black hole. The concept of BLR stratification arose from reverberation studies. Peterson & Ferland (1986) reported variability in the profiles of HeII4686Å and Hβ in NGC 5548, most notably the appearance of broader line components while the continuum luminosity increased.
Measuring the asymmetry of broad emission lines and including the indices in a statistical analysis therefore offers an approach to investigate how the relative kinematics of as-sumed layers in the BLR impacts the diversity of observed Type 1 AGN. One can also use the asymmetry index to trace the relative displacement of narrow emission lines and their shifted wings. This is particularly useful to trace the presence of outflows in the NLR (e.g., Heckman et al. 1981;Zamanov et al. 2002;Zakamska et al. 2016;Wang et al. 2018;Rakshit & Woo 2018 ). In this case the displacement of peak and base components is measured for the full profile (i.e. narrow core and shifted wings).
The asymmetry index was not included in C19 and was determined here starting from the reconstruction of the lines, using the Gaussian models listed in the catalogue. In this work we derived the following asymmetry indices: • ∆λ Hβ for the broad component of Hβ. Following C19, we retain the Gaussians with FWHM > 800 km s −1 as broad components. From the four Gaussian components used to fit the full Hβ profile, one was fixed by C19 at FWHM < 800 km s −1 to force the fit of the narrow emission line. Up to three Gaussians were thus used to fit the broad profile.
• ∆λ [OIII] for the full profile of [OIII]λ5007Å. Here both narrow and broad components are used for the measure.
We use the same method to measure these two indices. The description of ∆λ for any emission line presented here was first introduced by Heckman et al. (1981), and is further detailed in Winkler & Chauke (2014).
This method requires the identification of the wavelengths at which the line reaches 15% and 80% of its maximum intensity, which required the reconstruction of the multi-Gaussian line fits as described in Appendix A. We determined the maximum of the full line profile (or alternatively the broad components for Hβ) and measure 15% and 80% of this maximum. A schematic view of this measurement is presented in Fig. 4. The asymmetry index (∆λ) is defined by: Here λ b,high and λ r,high are measured at 80% fractional intensity and λ b,low and λ r,low at 15% fractional intensity. λ c,high is the wavelength of the center of the line at 80% fractional intensity, i.e., the midpoint between λ b,high and λ r,high . The parameters α b and α r are the distances from λ c,high to the blue and red edge of the line at 15% fractional intensity.
A negative (positive) ∆λ corresponds to a blue-ward (red-ward) asymmetry, and its value ranges from -1 to 1. This measure is not sensitive to the continuum emission (cf. correlation matrix shown in Section 4.1). We propagated the uncertainties of the FWHM of the fits into the measure of asymmetry and obtain typical errors of σ ∆λ ∼ 0.05. An important caveat is that using the uncertainties of the FWHM from the full profile, does not account for the lower, distinct and possibly more contaminated line components.
Examples of fit reconstructions of red-and blueasymmetric Hβ from our sample are displayed in Appendix B, alongside their Gaussian decomposition.

IMPACT OF FIT CONTAMINATION ON THE Hβ AND [OIII] ASYMMETRY INDEX DISTRIBUTION
The distributions for Hβ and [OIII] asymmetry indices are shown respectively in Fig. 5-6. The indices clearly follow very different distributions for the two emission lines. The following statistics were derived for the distributions of asymmetries with base centroids measured at 25% fractional intensity (i.e., the blue histograms in Fig. 5-6). The asymmetry indices of the [OIII] are predominantly blue-ward (skewness of distribution: γ = −0.54), which is consistent with the presence of a blue-shifted wing detected in a two-Gaussian model. The mean value of ∆λ [OIII] is ∼ −0.10. 1546 sources with [OIII] asymmetry below zero and 547 above zero are counted, confirming the excess of blue asymmetric [OIII] lines. The distribution of the Hβ asymmetry appears positively asymmetric (skewness of distribution γ = 0.36). The mean value of ∆λ Hβ is +0.03. An excess of red-asymmetric Hβ emitters is clearly measured (1244 sources out of 2100 with ∆λ Hβ > 0). The fractional intensity at which the lower velocity shift is measured for the asymmetry index has important repercussions on the overall distribution of asymmetries as Fig.  5-6 demonstrate. The distribution of ∆λ Hβ and ∆λ [OIII] are displayed for different fractional intensities at which λ b,low and λ r,low are set : 5%, 10%, 15%, 20%, 25%. While we always detect an excess in red-ward asymmetry for the Hβ profiles, the lower the percentage at which λ b,low and λ r,low are measured, the more the asymmetry index distribution is bimodal, with a clear separation between the symmetric and the red asymmetric Hβ modes. It is unclear if this red-ward excess in Hβ asymmetries is influenced by [OIII] and continuum contamination in the fit. The bimodality of the ∆λ Hβ distribution measured at lower base intensities, however, clearly indicates a preferred displacement between peak and base component which could arise from the near systematic contamination by [OIII] or its blue wing.
A similar effect is observed for the distribution of the [OIII] asymmetries: it becomes double peaked, with a distinct second mode appearing blue-ward of the symmetric sources for lower base intensities. For the measurement of the asymmetry index, we set the base of the lines at 15% fractional intensity, in order to avoid the detection of overlaying line models. This choice is motivated by two aspects of the distributions seen in Fig. 5-6 : 15% fractional intensity, the peak of the symmetric sources is more sensitive to low intensity, shifted profile components are detected by the measurement than at 25% and 20%. Furthermore, for this measurement configuration (peak at 80% and base at 15% fractional intensity), the bimodality which we identify as a signature of model degeneracy is not yet detected.
To further investigate the risk of model degeneracy, the dependency of the [OIII] asymmetry index on the shifts of the Hβ base component c 15 (the centroid shift at 15% fractional intensity with respect to the rest-frame wavelength) is presented in Fig. 7. We clearly observe an excess of sources in the quadrant of blue-asymmetric [OIII] lines and redshifted Hβ base components. The coloured markers correspond to high signal-to-noise ratio sources, which also appear to preferentially populate the lower right quadrant in the figure. 1196 sources out of 2100 (57%) have ∆λ [OIII] < 0 and c 15 > 0.
In the high signal-to-noise regime, i.e. S/N > 35, 278 out of 498 sources (56%) have ∆λ [OIII] < 0 and c 15 > 0. A binomial test showed that the clustering in the lower right quadrant is indeed significant. For both, c 15 and ∆λ [OIII] , the null hypothesis that the sources distribute equally on both sides of zero was rejected with p-values p < 10 −67 .
This result underlines the possibility of model degeneracy in the Hβ region. In their study of Extremely Reddened Quasars (ERQs), Perrotta et al. (2019), recently presented rare sources for which Hβ and [OIII] appear to blend (see their Fig. 1). Disentangling the presence of a shifted Very Broad Component from fit contamination would offer a cleaner window on the kinematics of the BLR, however, the fitting procedure in C19 did not allow for such a decomposition.

STATISTICAL ANALYSIS
This section presents the statistical analysis performed on the selected optical properties of our sample of 2100 SPI-DERS AGN, listed in Section 2.

Direct correlation
Our initial step is to generate a Spearman rank correlation matrix for our sample; the resulting matrix is shown in Fig. 8. Statistically insignificant correlations, defined by imposing a threshold of p max = 0.05/55 on the p-values, were masked. The threshold p max applies the Bonferroni correction for multiple statistical hypothesis testing (e.g., Haynes 2013).
7 Throughout this work, ρ denotes correlation coefficients, while p stands for p-value   The filled circles mark the sources with S/N above 35. The green dotted horizontal lines indicate typical 1σ uncertainties in the asymmetry parameter. The density contours containing 68%, 86%, 95% and 98% of all the sources in our sample were obtained from a Gaussian kernel density estimate.
Positive (negative) coefficients indicate a positive (negative) monotonic relation between two parameters. We have included a dendrogram (tree diagram) which clusters our data hierarchically using correlation as distance metric (for a review see Baron 2019). The clustered correlation matrix reveals structure in our parameter-subset.
• Similarly, an even stronger anti-correlation (ρ = −0.61) is measured between the equivalent width of [OIII] and r FeII , in agreement with the main EV1 anti-correlation.
• The X-ray and continuum luminosities, L X and L 5100Å , are related to the emission line properties in a similar manner. There is a strong correlation (ρ = 0.54 for L 5100Å and ρ = 0.55 for L X ) with the equivalent width of Hβ. The correlation between L 5100Å and W(Hβ) is consistent with photoionisation models in which increasing emission of the central engine results in a more luminous BLR. L X and L 5100Å decrease strongly (ρ = −0.44 and ρ = −0.51 respectively) with increasing flux ratio F(OIII)/F(Hβ).
In order to study the correlations with the ∆λ Hβ parameter more closely we perform a partial correlation analysis of FWHM Hβ , W([OIII]), r FeII and ∆λ Hβ (e.g., Baba et al. 2004). For this exercise, we first generate the Pearson correlation coefficients for this subset of parameters, then measure the strength of each of these correlations while controlling for the other confounding variables. In the simplest case, when only one confounding variable is accounted for, this is achieved by performing a linear regression for each of the two test variables with the confounding variable. The partial correlation of the two test variables is then obtained by measuring the correlation of the residuals resulting from the linear regressions. The results are shown in Fig. 9. When we control for FWHM Hβ and r FeII , the correlation between W([OIII]) and ∆λ Hβ significantly decreases (we measure a drop ρ P = 0.20, p < 10 −7 to ρ P, partial = 0.017, p = 0.45). The p-value of the latter partial correlation is above the corrected significance threshold: p = 0.05/6. Similarly, when we control for W([OIII]) and r FeII , the strength of the FWHM Hβ vs. ∆λ Hβ correlation decreases (ρ P = 0.28, p < 10 −6 to ρ P, partial = 0.075, p ∼ 10 −4 ). The anti-correlation between ∆λ Hβ and r FeII is, however, less affected when we control for FWHM Hβ and W([OIII]) ( ρ P, partial = −0.44, p < 10 −6 instead of ρ P = −0.52, p < 10 −6 ). These results demonstrate that ∆λ Hβ parameter is principally related to the relative strength of the iron emission and is marginally linearly correlated to FWHM Hβ , indicating that the broader profiles tend to be more red-skewed, as also suggested by Marziani et al. (2013a).
The key result so far is not the rediscovery of correlations mapped out by e.g., Boroson & Green (1992); Grupe (2004); Shen & Ho (2014), but the fact that these relations hold for type 1 AGN in general, up to a redshift limit of at least z = 0.80 and for luminosities up to L X ∼ 10 46 erg.s −1 .
A new insight provided by the correlation matrix is the correlation behaviour of ∆λ [OIII] with parameters related to the BLR emission : • The asymmetry index of [OIII] appears to be marginally related to the rest of the parameters, suggesting the absence of kinematic linkage between the inner region of the AGN and the NLR. This result is in contrast to Zamanov et al. (2002), who reported evidence for correlation between the shifts of the high ionisation lines CIVλ1549Å and the shifts of the [OIII] lines, inferring a possible linkage between the NLR and BLR.

Principal components
The second step consisted of running a PCA on the standardized data set to determine which parameters are contributing most to the total variance of our sources. It is particularly important to scale the parameters to unit variance since PCA is sensitive to the variance of the parameter  distributions. PCA is an orthogonal linear transformation which is often used as a dimensionality-reduction algorithm (for a review see Jolliffe & Cadima 2016). It yields the eigenvectors in parameter-space which point in the direction of total maximal variance in the dataset. For an N × M dataset, the first component is found by minimising the distance between points in the M-dimensional parameter space and their orthogonal projections onto an M-dimensional vector, which simultaneously increases the variance of the projected Figure 10. The component coefficients (factor loadings) defining the first and second principal components : EV1 (green bars) and EV2 (red bars). EV1 appears to be heavily dominated by the anti-correlation of FWHM Hβ and the strength of the [OIII] emission with r FeII . A relatively large linear coefficient links the Hβ asymmetry to EV1. EV2 is heavily dominated in equal measure by the X-ray and optical luminosities.
points. The second component is chosen to be orthogonal to the first and is determined in the same manner given this condition. The associated eigenvalues measure the amount of explained variance in each principal component (eigenvector). The variables (our parameters) are linked to the principal components by linear coefficients. These coefficients yield the amount of a variable's variance explained by the component. 8 While PCA is a very common tool in astronomy it comes with several disadvantages, such as the linearity of the dimensionality reduction, which might not capture the full complexity of our data. One must in general be careful with the choice of the particular dimensionality reduction one uses and its interpretation (e.g., see Baron 2019).
The first two principal components are presented in Fig.  10. The bar diagrams show the values of the correlation coefficients which link the parameters (our variables) to the given component. The first and second principal component respectively explain 25.2% and 20.3% of the total variance in the selected parameter subset.
The first principal component (Eigenvector 1, EV1), is anti-correlated to the strength of the iron emission and correlated to FWHM Hβ and the (relative) strength of the [OIII] emission, i.e., the diversity in the selected optical features of our SPIDERS AGN sub-sample is dominated by the anti-correlation of W([OIII]) and W(FeII). The Balmer profile asymmetries ∆λ Hβ are noticeably correlated to EV1, which corroborates the findings of Z10, who stress that the asymmetry parameter could in essence be used as surrogate 4DE1 parameter. The authors note, however, that ∆λ Hβ is also substantially orienting Eigenvector 2 (EV2). This behaviour is not observed in our PCA results. Our EV1 is consistent with early results in Grupe (2004).
The second principal component, EV2, is strongly dominated by the X-ray and optical continuum luminosities. The flux ratio F([OIII])/F(Hβ) is correlated to EV2. The equivalent width of Hβ once again significantly contributes to the orientation of the principal component. EV2 is also consistent with Grupe (2004).

BLACK HOLE MASS AND EDDINGTON RATIO
We extend the analysis to derived AGN properties: the black hole mass and the Eddington ratio. In Boroson & Green (1992) EV1 was related to the Eddington ratio L/L Edd , tracing the accretion power of the observed SMBH. Their EV2 was largely dominated by luminosity and related to the Baldwin effect (Baldwin 1977;Baldwin et al. 1978;Dietrich et al. 2002). The Eddington ratio has indeed long been considered a prime candidate to explain the observed diversity in optical AGN features (e.g., Sulentic  The FWHM of broad emission lines serves as proxy for the line of sight velocity of the BLR gas, which in combination with the distance to the black hole yields an estimate for its mass : f BLR denotes the geometric form factor of the BLR. The distance r BLR from the SMBH to the line-emitting region is obtained from the delay in response of the line fluxes to continuum changes. The FWHM BLR of a line from the emitting region corresponds to the radial velocity component of the BLR gas. The Keplerian velocity v Kepl of the gas is obtained by correcting for projection effects encoded in the geometric form factor : v Kepl 2 = f BLR FWHM BLR 2 . One can construct an expression of f BLR , which takes in account the source orientation θ and the shape of the velocity field described by the motion of the gas in the BLR (e.g., equation 12 in Negrete et al. 2018), by making assumptions about the structure of the disk.
The black hole masses and the Eddington ratio for our sample are derived from Hβ using the calibration developed in Assef et al. (2011), who based their estimates on the BLR radius-luminosity relation of Bentz et al. (2009). L λ corresponds to the monochromatic luminosity at 5100Å and the parameters A, B and C are derived from reverberation mapping studies: A = 0.895, B = 0.52 and C = 2. From these black hole masses, C19 provided estimates of the Eddington ratio following: where G is the gravitational constant, c the speed of light, m p the proton mass and σ T the Thomson scattering cross-section. The bolometric luminosities L Bol used here are obtained from the bolometric corrections presented in Richards et al. (2006) : L Bol = 9.26L 5100Å . This factor is in principle dependent on L 5100Å , but for the purpose of this work, this assumption is not crucial. Fig. 11 displays how the Eddington ratio and black hole mass scale with the projection of the sources onto our newly determined EV1 and EV2.
We can strongly corroborate that EV1 correlates with the Eddington ratio. An even stronger anti-correlation is that of EV2 with black hole mass. However, the log M BH -EV2 relation breaks down at lower redshifts (z ∼ 0.2). In order to determine to which extend this correlation is a sideproduct of the L 5100Å and FWHM Hβ principal coefficients in EV2, we once again perform a partial correlation test. We start by measuring the Spearman rank correlation of EV2 with black hole mass, while controlling for L 5100Å and FWHM Hβ . In Fig. 12 the Spearman correlation coefficients measured for log M BH , EV2, L 5100Å and FWHM Hβ are displayed in the left panel. In the right panel of Fig. 12, the partial correlation coefficients results are presented. As already discussed in Section 4, this exercise consisted in performing a (multi) linear regression of EV2 and log M BH with the confounding variables. In this version we compute the Spearman correlation coefficients of the residuals of these two multilinear regressions. EV2 is still strongly anti-correlated to the black hole mass (ρ = −0.56). This is not surprising since in Eq. 4, log M BH is a non-linear function of FWHM Hβ and L 5100Å , while the Eigenvectors obtained from PCA are obtained through linear orthogonal transformations of the initial parameter. The monotonic anti-correlation measured in the residuals might thus simply arise from the non-linearity of Eq. 4. More generally, this result implies that one can construct a linear combination of the parameters which dominate the orientation of EV2 (see Fig. 10) in order to estimate the black hole mass at redshifts z > 0.2.
We propose that EV2, through its anti-correlation to black hole mass, might be related to the evolution of the broad line AGN population. Over the observed redshift range, the median black hole mass increases with decreasing redshift. This trend a combination of two effects: the 2RXS flux limit resulting in the Malmquist bias and the decreasing number density of high black hole masses at lower redshifts. This downsizing of black hole masses across cosmic time remains a matter of debate (cf. downsizing models, e.g., Fanidakis et al. 2012).
We conclude that the Eddington ratio and the black hole mass are related to the principal components which explain ∼ 45% of the total variance in our data.

ASYMMETRY OF THE BROAD Hβ EMISSION LINE; A MARKER OF TYPE 1 AGN DIVERSITY
EV1's direction in parameter space has been used in the past to establish an analog to the stellar Hertzsprung-Russel diagram for the Type 1 AGN population. In the plane spanned by FWHM Hβ and r FeII (i.e., the first two dimensions of the 4DE1) the domain occupied by Type 1 AGN has been pre-  sented as a quasar main sequence (Marziani et al. 2001;Sulentic et al. 2011). We investigate the role played by Hβ line asymmetries in this framework. Since they appear to be one of the central sources of variance in our sample, they can be source of optical diversity of Type 1 AGN. Fig. 13 presents the sample in the FWHM Hβ vs. r FeII plane, also called EV1 plane. The sources have been separated in two subsets: blue and red asymmetric Hβ emission. More specifically these two sub-samples were constructed ac-cording to : ∆λ Hβ > 0.07 and ∆λ Hβ < −0.07. These criteria exclude symmetric sources and account for the typical uncertainties of the asymmetry index (σ ∆λ Hβ ∼ 0.05). Contours of the bivariate Gaussian kernel density estimates (KDE, e.g., Silverman 1986) for each of the subsets are also indicated. The contours delimit the areas containing 68%, 86%, 95% and 98% of the data points. 9 The sample occupies an L-shaped form in the plane. Using simulations, C19 demonstrate that the absence of high FWHM Hβ and high r FeII sources is to an extent due to limitations of the spectral fitting method. Most of the sources with blue asymmetric Hβ emission profiles and those with red-asymmetric profiles appear to occupy different domains of this sequence. While the AGN with blue-ward asymmetric Hβ spread over the full r FeII range and a large portion of the FWHM Hβ range, the red-ward asymmetric sources are concentrated at lower values of r FeII , while dominating the higher segment of FWHM Hβ . This result is essentially suggesting that high accretors, with relatively large r FeII values (e.g., Panda et al. 2019), show mainly blue-ward asymmetric Hβ profiles. We also observe a considerable overlap of red and blue asymmetries for sources with moderate FWHM Hβ and r FeII .
In their in-depth study of Hβ line profiles of ∼ 470 low-z SDSS (DR5) quasars, Z10 show that low accretion rates sources (low r FeII ) possess a typically red-asymmetric Hβ profile, while high accretion rate sources tend to prefer blue-asymmetric Hβ profiles. Figs. 14-15 can directly be compared to Fig. 6a and 9a presented in their work: • The size of our sample allows us to extend their results by affirming that blue asymmetric Hβ profiles can be found Figure 13. The horizontal trend of the H β asymmetry in the EV1 plane is made clear in the figure above. The contours derived from the KDE, delimit the areas in which 68%, 86%, 95% and 98% of the data points are confined. These contours reveal two very different occupation domains for red-and blue asymmetric Balmer emitters : sources with red-asymmetric Hβ are confined at low r FeII values, while sources with blue asymmetric Hβ seem to extend over the full EV1 sequence.
in sources with relatively high FWHM Hβ (∼ 8000 km s −1 ) and along the full EV1 sequence. There is, however, a clear decrease in blue asymmetric Balmer profiles at FWHM Hβ ∼ 4000 km s −1 , consistent with the low-redshift separation in Population A/B introduced by Sulentic et al. (2000a). An excess of red-asymmetric Hβ profiles is indeed observed for larger widths. Sources with lower FWHM Hβ do not, however, display more symmetric profiles. We thus argue that there is no evidence for a systematic relationship between FWHM Hβ and ∆λ. The excess of red asymmetric Hβ emission for larger widths might be the signature of an additional redshifted broad emission component.
• Fig. 15 strongly confirms the decreasing trend presented in figure 5 of Boroson & Green (1992) and figure 9a of Z10. Combined with the Spearman correlation coefficient of ρ ∼ −0.5 reported in the previous section, we can state that, up to redshift of z ∼ 0.8, low accretors tend to show more red asymmetric profiles, with a significant decrease of red-asymmetric Hβ emitters beyond r FeII ∼ 1.5. While blueasymmetric profiles are reported over the full r FeII range, the shifts become more pronounced at higher values of the iron-Hβ flux ratio.
Z10 discuss the VBLR-emission as potential origin of typically red asymmetries in Hβ (2.2). In a layered model of the BLR, the asymmetry index traces the displacement between the the Hβ broad component and its very broad component. Hβ lines often show a broad, redshifted, lowintensity feature in their profiles (Peterson & Ferland 1986;Corbin 1995;Brotherton 1996). The presence of a distinct emission region, the VBLR (Popović et al. 2004;Marziani et al. 2010), has been proposed to explain the characteristic red wing of broad Balmer emission. The VBLR is expected to be located in the inner region of the AGN: The typical  widths of the very broad component in Hβ indicate much higher Keplerian velocities than the classical broad line region. The presence of a VBC in the profile shapes has to be corrected for when broad low ionization lines are used as virial broadening estimator (Marziani et al. 2013a). In Appendix B, three sources of our sample, which have one broad Gaussian component with FWHM > 10.000km s −1 , are presented.
In order to further characterise the kinematics of two potentially distinct emitting regions, we compute the centroid shifts of the Hβ profiles at 15% (c 15 ) and 80% (c 80 ) fractional intensity. Fig. 16 displays the correlations between these shifts and the asymmetry index of Hβ, colour-coded according to the FWHM Hβ . The shifts with respect to the rest-frame at 80% fractional intensity are negatively related to the asymmetry index, i.e., the more blue shifted the top of the profile, the more red-ward asymmetric the complete profile and vice-versa. Symmetric profiles range over the full shift range.
For the centroid shifts at c 15 fractional intensity, an opposite and stronger trend is observed. c 15 correlates positively with ∆λ Hβ : a red-ward asymmetric Hβ profile can be associated to a redshift of the broad component. As expected, sources with symmetric Hβ span over the full c 15 range. The highest values of FWHM Hβ occur for red asymmetric Hβ consistent with the results in Fig. 14. A clear picture arises from the c 15 distribution in the EV1 plane (Fig. 17). The absolute value of the shifts at the profile base appears to decrease along the EV1-sequence, while c 15 shows clear trends with both r FeII and FWHM Hβ . This picture is consistent with the presence of a redshifted VBLR.
The base shifts of the Hβ profiles in our sample are tightly correlated to the asymmetry index ∆λ Hβ . The largest shifts are observed for the most red-asymmetric profiles. The highest base component redshifts are found at the top of the Type 1 AGN main sequence, i.e., at the highest widths of Hβ: the broadest components in our sources' Hβ profiles are preferentially redshifted. The correlation of c 15 with EV1 provides additional information: as can be seen in Fig. 18, the positive base centroid shifts appear to correlate with EV1. Interestingly, the broad redshifted Hβ components are found for the sources with the highest black hole masses.
The broad component redshifts are thus a source of variance in our sample. If the redshifted, very broad components are interpreted as the signature of a VBLR, we can argue that the presence or absence of a distinct, inner emitting shell in the BLR (or more precisely its redshift) is a strong source of diversity in our sample of Type 1 AGN. Furthermore we report that our results are consistent with Marziani et al. (2009), who find a systematic increase of FWHM Hβ with source luminosity. More precisely, Marziani et al. (2009) find a scaling relation between the luminosity and the contribution of the VBLR to the full Hβ flux in population B (> 4000 km s −1 ). The absence of red asymmetric Hβ lines at higher r FeII shown in Fig. 15 might be due to strong accretion disk winds preventing the formation of a VBLR.
In Appendix C, we explore a toy model of BLR obscuration, which accounts for the presence of a VBLR. If we consider that blue-and red-ward asymmetric Hβ profiles are signatures of two distinct kinematic states of the BLR, we should divide our sample according to the sources' asymmetry indices to investigate the properties of each subsample. As in the previous section, we formed two subsamples with ∆λ Hβ < −0.07 (blue asymmetric) and ∆λ Hβ > Figure 16. The Hβ centroid shifts at 15% and 80% fractional intensity as a function of Hβ asymmetry. c 15 appears more tightly related to ∆λ Hβ than c 80 . The kernel density contours show the sub-sample of sources, for which Hβ has not been fit with broad Gaussians (FWHM < 10.000 km s −1 ). Typical 1σ uncertainties in the measurement of ∆λ Hβ are indicated as vertical red dotted lines.
0.07 (red asymmetric). We constructed a Spearman rank correlation matrix for the parameter subset, as defined in Sec. 2.1, of each sub-sample and directly compared their correlation spaces. A striking difference is observed for the correlation of the optical and X-ray luminosity and the iron emission. For blue asymmetric Hβ emitters, L X and L 5100Å correlate positively with the equivalent width of the iron emission W(FeII) (Spearman correlation coefficients: ρ S = 0.47 and ρ S = 0.55, respectively).The sources with red asymmetric Hβ have a much weaker correlation between the source luminosities and the equivalent width W(FeII) (ρ S = 0.023 for L X and ρ S = 0.097 for L 5100Å ). The correlation of L 5100Å and W(FeII) for red-asymmetric sources is p = 0.53 and can safely be considered as insignificant. This contrasting behaviour is presented in the left panel of Fig. 19. In order to confirm the different clustering of red-and blue-asymmetric Hβ emitters in the plane spanned by L 5100Å and W(FeII), we  binned log W(FeII) in the range 1.0 − 1.9 and bootstrapped the sub-samples in each bin (10.000 resamples). We obtained the mean L 5100Å for each bin from the sampling distribution of the means. The 68% confidence intervals were derived using the percentile method (from the 16.0% and 84.0% percentiles). This method should yield relatively good estimates of the error, given the symmetric shape of the sampling distribution of the mean in each bin (shown in the right panel of Fig. 19). The positive correlation between L 5100Å and W(FeII) for sources with blue asymmetric Hβ is confirmed.
No trend is found for sources with red-asymmetric Hβ. The error bars, showing the standard deviation of the the bootstrap samples' means, do not overlap. The correlation of L 5100Å and W(FeII) for the blue asymmetric Hβ population withstood a partial correlation test, where we controlled for the confounding W(Hβ). Similarly, when we control for the FWHM FeII of the Gaussian kernel used to fit FeII, L 5100Å and W(FeII) remain positively correlated. The dependency of the FeII flux on the continuum luminosity in blue-asymmetric Hβ emitters is only marginally related to the broadening of the iron lines. This contrasting behaviour also manifests itself in the correlation of the luminosity parameters with the flux ratio r FeII . Blue-asymmetric asymmetric Hβ sources are significantly correlated r FeII and source luminosities (ρ S = 0.20 for L 5100Å and ρ S = 0.19 for L X ). For red asymmetric Hβ sources, these quantities are anti-correlated (ρ S = −0.16 for L 5100Å and ρ S = −0.25 for L X ).
Another noticeable difference is observed for the ∆λ [OIII] parameter. In the sub-sample of blue-asymmetric Balmer emitters, ∆λ [OIII] and F([OIII])/F(Hβ) are anti-correlated (ρ S = −0.26). If blue shifted wings are interpreted as the signature of NLR outflows, this result would imply that outflow velocities increase with increasing [OIII]/Hβ flux ratios. In the red-asymmetric sub-sample, these two parameters are positively correlated (ρ S = 0.17).
Blue-ward asymmetries of low-ionisation lines have been related to radiation driven outflows (Marziani et al. 2013b). They occur in the high r FeII bins of population A. (e.g. Ganci et al. 2019), which contain the highest accretion rate sources along the main sequence (e.g. Sun & Shen 2015;Sulentic et al. 2017;Panda et al. 2019). The origin of FeII emission in the BLR has been a long-standing matter of research. Several lines of evidence support an emission of FeII in the outer parts of the BLR, while Hβ might be emitted closer the the black hole (e.g. Rodríguez-Ardila et al. 2002;Barth et al. 2013;Marinello et al. 2016). The physical conditions for iron ionisation have been investigated in detail. Wills et al. (1985) showed that photoionisation models such as the locally optimally emitting clouds (Baldwin et al. 1995) might not suffice to account for the total observed FeII flux in AGN. For a review of necessary conditions for the formation of the observed FeII in photoionisation models see Collin & Joly (2000). Models of collisional excitation (e.g. Baldwin et al. 2004;Joly et al. 2007), in addition to continuum and line fluorescence (e.g. Sigut & Pradhan 1998;Marinello et al. 2016), have been considered. These models make predictions on the physical conditions of the FeII emitting region, such as shielding from the continuum source and high densities.
In a flattened, horizontally stratified cloud distribution in Keplerian motion around the central black hole, the difference in sensitivity to continuum of the FeII emitting region might arise from different degrees of exposure to the central continuum source of the dense clouds from which the iron is emitted. Radiation driven winds produce outflows of Hβ emitting gas, thereby exposing the previously shielded FeII regions more directly to continuum emission. The scaling between W(FeII) and the L 5100Å might thus be due to the contribution of UV fluorescence to the excitation of Fe+ levels. Investigating FeII fluorescence in HII regions, Rodríguez (1999) uses the sensitivity to fluorescence of [FeII]λ4287, as well as the continuum insensitive [FeII]λ8617 to investigate the role played by the UV radiation field in the formation of this line. If the presented model of unshielding holds, the intensity ratio I([FeII]λ4287)/I([FeII]λ8617) should scale differently with the UV/optical continuum luminosity for the red-and blue-asymmetric Hβ populations. This may be investigated in further work.
In this scenario, the distance of the FeII region to the SMBH should not depend on the presence of outflows. We performed a two-sample Anderson-Darling test on the distributions of the FWHM FeII of the Gaussian kernel which was convolved with the FeII template in the red-and blueasymmetric Hβ subsets. 55 sources in our sample have no FeII flux detection and were not considered for this test. The obtained p-value is p > 0.25. We thus cannot reject that the FWHM FeII for red-and blue-asymmetric Hβ sources were sampled from the same distribution. Under the assumption of a Keplerian velocity field, the radial distribution of the FeII emitting clouds is not affected by the presence of outflows.
We conclude that our data is consistent with a flat, stratified and self-shielding model of the BLR, in which the FeII and Hβ emissions originate at different radii. When the inner clouds are driven out of the plane by outflows during strong accretion events, the previously shielded, neutral and dense FeII clouds are more exposed to the ionizing continuum radiation, which is clearly seen in the FeII fluxluminosity scaling.

Evidence for model degeneracy: FeII vs. Hβ
We note that C19 fitted the emission of the iron complex over the full Hβ centered region, using the I Zw 1 template (Boroson & Green 1992). Since all the emission features of FeII in the spectral window were taken in account during the fits, the effect of contamination by Hβ should be marginal on the total FeII flux, i.e. W(FeII) should not be strongly affected by increasing Hβ line flux. However we can investigate if the FeII contaminates the Hβ fit on its blue side.
For blue asymmetric Hβ emitters, we seek to test if the the equivalent width of FeII correlates with the Hβ asymmetry index, i.e. we want to determine if the excess flux on the blue side of Hβ is due to contamination by iron.
On the sub-sample with, ∆λ Hβ < −0.07, we measure the Spearman rank correlation coefficient of | ∆λHβ | and W(FeII). We obtain ρ = 0.087 with a p-value of p = 0.08. Even if we do not correct the significance threshold for multihypothesis testing, one can state that there is no significant correlation between ∆λ Hβ and W(FeII). We note however that the equivalent width of Hβ is a confounding variable for both: ∆λHβ (p = −0.26) and W(FeII) (p = 0.29). We thus have to once again perform a partial correlation analysis, marginalising over W(Hβ) in order to obtain the real correlation behaviour of the blue-asymmetries and the iron strength. The partial correlation coefficient is p P = 0.17 with a p-value of p ∼ 10 −4 . There is a significant correlation between ∆λ Hβ and W(FeII) once we control for W(Hβ). We thus cannot exclude that the iron complex contaminates the blue wing of ∆Hβ.

CONCLUSIONS
In this work we have performed a statistical analysis of a high signal-to-noise ratio (S/N > 20) sub-sample of the SDSS-IV/SPIDERS DR14 VAC of X-ray selected Type 1 AGN. The sample included 2043 sources spanning an X-ray luminosity range of L X = 1.9 × 10 41 − 9.9 × 10 45 erg s −1 , up to redshift z = 0.80. It is the largest sample of high signal-tonoise ratio sources used for the PCA-based, statistical study of Type 1 AGN spectral properties to date.
• We used PCA as a central tool to determine the source of variance in our data and have mapped a correlation space which is remarkably consistent with previous studies of the optical Eigenvector 1 (e.g., Boroson & Green 1992;Sulentic et al. 2000a;Grupe 2004;Shen & Ho 2014), while probing a larger cosmological volume.
• We confirm that the Eddington ratio and the black hole mass are significantly related to the observed diversity of Type 1 AGN, through their correlation to EV1 and EV2.
• Studying Hβ line shapes in this context, we find blue asymmetric emission profiles for the full Type 1 AGN main sequence, while red asymmetries only appear for low accretors. Z10, investigating Hβ profile asymmetries in the 4DE1 context, have suggested that profiles with lower FWHM Hβ tend to be more symmetric while profiles with larger FWHM Hβ are preferentially red-asymmetric. The larger number of sources in our sample enabled us to complete this picture: While we do observe a larger portion of sources with red asymmetry index at high FWHM Hβ , lower width profiles appear to cover the full range of asymmetries observed in our sample, i.e., lower FWHM Hβ do not lead to more symmetric profiles. We can, however, confirm a strong trend between r FeII and ∆λ Hβ . In their discussion of the physical origins of Hβ profile shapes, Z10 identify VBLR-emission and disk winds as main candidates.
• A sub-class of our sources does indeed show very shifted, broad components in their Hβ emission, possibly due to the presence of a distinct emitting gas distribution in the inner regions of the BLR. The redshift of this VBLR correlates with Eigenvector 1 and might thus be related to Type 1 AGN diversity. • Exploring parameter correlations for blue-and redward asymmetric Hβ emitters separately, we observed that FeII line flux correlates differently with source luminosity for redand blue-ward asymmetric Hβ emitters. We discussed this effect in the light of a flattened, self-shielding BLR, in which the FeII emitting clouds are located are larger radii than the ones emitting the Balmer lines. The outflows, we associate Hβ blue-ward asymmetries with, might in such a configuration strip the FeII region of its shield, i.e. the innermost parts of the BLR. However, we find evidence for inter-line contamination between FeII and Hβ, which might play a confounding role in this effect.  APPENDIX C: TESTING A SIMPLE OBSCURATION SCENARIO Richards et al. (2002) proposed that the spatial configuration resulting from the BLR's inclination and obscuration might explain the asymmetries in the observed profile of CIVλ1549Å (see also 4.3. in Z10). In this picture, the asymmetry arises from the attenuation of redshifted photons by obscuring material. Similarly Gaskell & Harrington (2018) developed a model of outflowing clouds of dust, blocking the line of sight to the inner disk-shaped BLR, which predicts broad emission profiles similar to those observed. An alternative view of the effect of the BLR and outflows is developed by Czerny et al. (2017) (see also Czerny & Hryniewicz 2010). In this model the dusty clouds are driven out of an optically thick disk, and are eventually exposed to radiation by the central source. The dust evaporates, removing the matter/radiation interaction and results in the fallback of gas clouds. This scenario is referred to as Failed Radiatively Accelerated Dusty Outflow (FRADO). The physics of dust sublimation and the consequences for the BLR disk structure were further investigated by Baskin & Laor (2018). Fig. 16 allows us to investigate a simple obscuration scenario in which the observed asymmetry of the line shapes of Hβ would result from partial obscuration of a flattened BLR in Keplerian motion around the black hole. If the receding or incoming region (with respect to our line of sight) of a disk like structure is obscured by optically thick material, the emitted broad lines would be respectively red-or blue asymmetric (see for example Fig. 2 in Gaskell & Harrington 2018), since one peak of the double-peaked emission profile would be attenuated. In this scenario, the low-intensity, 'excess' emission measured on blue (or red) side of the line would correspond to the attenuated blue (or red) horn of the disk-like emission.
Furthermore, if we assume no further kinematics in our BLR model, the attenuation of the red of blue peak of the disk like emission would lead to a centroid shift configuration : • Blue-ward asymmetric Hβ (∆λ Hβ < 0): Redshift of the peak component (c 80 > 0) and blueshift of the base component (c 15 < 0) • Red-ward asymmetric Hβ (∆λ Hβ > 0): Blueshift of the peak component (c 80 < 0) and redshift of the base component (c 15 > 0) Fig.16 supports this scenario. The results in the previous section suggest that the broad Balmer shapes might be affected by the presence of a distinct emitting region: the VBLR. The next step is to test if the obscuration scenario is compatible with the presence of a very broad component related to the stratification of the BLR. To this effect we identify sources for which the broad Hβ has been fit broad Gaussians. Fig. C1 presents the peak wavelength of all the Gaus-sians which have been used to fit the components of Hβ as a function of their FWHM. We clearly identify a sequence of very broad Gaussians at velocities > 10.000 km s −1 which are mostly redshifted with respect to systemic redshift. We investigated the reduced χ 2 distribution of the fits of the sources which have at least one Gaussian component of FWHM > 10.000 km.s −1 . The sources which have a very broad component are not biased towards higher reduced χ 2 (and thus higher BIC) values, indicating that they are not preferentially found in less secure fits. The statistical significance of a second mode in the distribution of the Gaussian FWHM in our sample at FWHM > 10.000 km s −1 was assessed with a Silverman test of multimodality (Silverman 1981, for its calibration see e.g. Hall & York 2001;Ameijeiras-Alonso et al. 2016). We could confirm that this second mode of very broad Gaussian components is indeed significant. The test is presented in details in Appendix D.
We empirically identify this sequence as the sub-class of sources which contain a VBLR. As expected the broadest Gaussians have the largest contributions to the full line luminosity. Fig.16 includes the density contours containing the sources for which Hβ has not been fitted with a very broad Gaussian (FWHM < 10.000 km s −1 , i.e., the sources to the left of the vertical dashed line in Fig.C1). The trend observed for all objects in our sample is preserved for this subset, supporting the possibility of a partially obscurated flattened BLR + VBLR model. The red-and blueshifts of the Hβ base and peak components are symmetrically distributed, which is consistent with the obscuration scenario. In a simple model in which wing attenuation via partial obscuration is the sole source of asymmetry of the broad emission lines, we expect the Hβ line fluxes to be symmetrically distributed for ∆λ Hβ < 0 and ∆λ Hβ > 0, as no wing of the emission line profile should be preferentially attenuated. The left panel in Fig.C2 shows the equivalent width of Hβ plotted against its asymmetry index for our full sample. We observe an offset in W(Hβ) between the population of red-and blue-asymmetric emitters. We argue that this offset might once again be the signature of the presence of a VBLR in some of our sources. If we display the kernel density contours for the subset of sources which were not fitted with a very broad Gaussian, a more symmetric distribution of equivalent widths is produced (i.e., line fluxes), as shown in the right panel of Fig. C2.
These simple tests have demonstrated that our sample is consistent with a model of partial obscuration coupled to a stratified BLR (one that might or might not contain a VBLR).

APPENDIX D: VBC MODE SIGNIFICANCE WITH A SILVERMAN TEST OF MULTI-MODALITY
The FWHM of broad Gaussian components (> 800 km s −1 ) that were fit to the Hβ emission lines show a bimodal distribution. The two modes appear to distinguish broad components (800 km/s < FWHM < 10.000 km/s) and very broad components (FWHM > 10.000km s −1 ). Out of the 8400 Gaussian components used to fit the 2100 sources of the sample, 4936 (58%) can be considered broad, with FWHM > 800 km s and 1080 (13%) are very broad components with FWHM >  10.000 km s. The significance of the second mode was assessed using the Silverman test of multimodality (Silverman 1981). We present the outline of this test here.
• We extract all the Gaussians used in the fits of our sample with the condition: FWHM > 800 km s −1 & λ peak 4830Å & λ peak 4900Å. This cut removes narrow Gaussians, as well as Gaussian which central wavelength was shifted to the edge of the fitting window. 4275 Gaussian components remain. We show the resulting density histogram of FWHM in Fig. D1.
• For a range of bandwidths h = [0.01, 2.00] (step size s= 100), we performed a Gaussian kernel density estimate (KDE) of the distribution of FWHM (centered at zero and scaled to unit variance). An example probability density function (pdf) from KDE is displayed in Fig. D1 (red line). For n observations in our sample, the bandwidth h smooths the Gaussian kernel estimate as follows: 2h 2 (D1) • We then detect the critical bandwidths h DATA c, j , at which new local maxima appear in the pdf, when h is decreased. For 1 < j < n (where n corresponds to the number of observations), the critical bandwidths h DATA c, j are formally defined as the minimum h, at which the Gaussian density estimate has no more then j maxima. In Fig. D2, the bandwidth h and the corresponding number of detected local maxima j are displayed. Critical widths are symbolized by vertical lines.
• In order to assess the significance of a mode detection, we follow Fischer et al. (1994) and set up a null hypothesis test, which is based on the observation that the number of modes off (x) decreases with increasing bandwidth h. If k is the true number of modes: Null hypothesis H 0 : k = j Alternative H 1 : k ≥ j + 1 for j = 1, 2...
For each h DATA c, j detected in our data, we perform a Gaussian kernel density estimate on 100 smoothed bootstrapped samples, using h DATA c, j as bandwidth. These resamples are obtained, by sampling smoothed FWHM s , using bootstrapping (with replacement). The exact computation of FWHM s is as follows: were FWHM r are resampled FWHM, σ, the standard deviation of the bootstrapped sample and is a random standard variable. For the null-hypothesis test Silverman (1981) motivates a p-value defined for a given number of maxima j as: , and thus each mode detection, we compute a p-value. We assess how many times the density estimates with fixed bandwidths of the 100 resamples have at most j detected modes. The ratio of these realizations and the total number of simulations are taken as p-value. The resulting p-values are listed in Table D1.
• Following Silverman (1981), we read this table as hierarchical successive significance test. As soon as for a new j, the p-values drop below an α significance threshold, we interpret the detection of j modes as significant. From Table  D1, if we adopt a significance threshold α = 0.05, it appears that the bimodality of the distribution of our Gaussian components FWHM is significant, i.e. the very broad mode of Gaussian components is significant. Figure C2. Left panel: The equivalent width of the broad Hβ emission as a function of its asymmetry Index. The blue and red density contours are included to improve the visualisation of the red-and blue asymmetric sub-samples. An offset between the two populations is observed. The horizontal red (blue) line mark the mean of the bootstrapped mean W(Hβ) of the the red-(blue-)ward asymmetric sources. The red and blue, horizontal dashed lines correspond to the 3σ confidence contours, as derived from the percentile method at 0.13% and 99.87% of the sampling distribution of the mean. Right panel: Same format as on the left. The blue and purple contours respectively represent the sources with blue-and red-asymmetric Hβ, for which the Gaussian models had FWHM < 10.000 km s −1 .
[h] Figure D1. Density of standardadised Gaussian component FWHM. This paper has been typeset from a T E X/L A T E X file prepared by the author. [h]