Extreme quasars as distance indicators in cosmology

Quasars accreting matter at very high rates (known as extreme Population A [xA] quasars, possibly associated with super-Eddington accreting massive black holes) may provide a new class of distance indicators covering cosmic epochs from present day up to less than 1 Gyr from the Big Bang. At a more fundamental level, xA quasars are of special interest in studies of the physics of AGNs and host galaxy evolution. However, their observational properties are largely unknown. xA quasars can be identified in relatively large numbers from major optical surveys over a broad range of redshifts, thanks to selection criteria defined from the systematic changes along the quasars main sequence. It has been possible to build a sample of about 250 quasars at low and intermediate redshift, and larger samples can be easily selected from the latest data releases of the Sloan Digital Sky Survey. A large sample can clarify the main properties of xA quasars which are expected - unlike the general population of quasars - to radiate at an extreme, well defined Eddington ratio with small scatter. As a result of the small scatter in Eddington ratio shown by xA quasars, we propose a method to derive the main cosmological parameters based on redshift-independent"virial luminosity"estimates from measurements of emission line widths, roughly equivalent to the luminosity estimates based from line width in early and late type galaxies. A major issue related to the cosmological application of the xA quasar luminosity estimates from line widths is the identification of proper emission lines whose broadening is predominantly virial over a wide range of redshift and luminosity. We report on preliminary developments using the AlIII 1860 intermediate ionization line and the Hydrogen Balmer line H-beta as virial broadening estimators, and we briefly discuss the perspective of the method based on xA quasars.

Quasars accreting matter at very high rates (known as extreme Population A [xA] quasars, possibly associated with super-Eddington accreting massive black holes) may provide a new class of distance indicators covering cosmic epochs from present day up to less than 1 Gyr from the Big Bang. At a more fundamental level, xA quasars are of special interest in studies of the physics of AGNs and host galaxy evolution. However, their observational properties are largely unknown. xA quasars can be identified in relatively large numbers from major optical surveys over a broad range of redshifts, and efficiently separated from other type-1 quasars thanks to selection criteria defined from the systematically-changing properties along the quasars main sequence. It has been possible to build a sample of ∼250 quasars at low and intermediate redshift, and larger samples can be easily selected from the latest data releases of the Sloan Digital Sky Survey. A large sample can clarify the main properties of xA quasars which are expected-unlike the general population of quasars-to radiate at an extreme, well-defined Eddington ratio with small scatter. As a result of the small scatter in Eddington ratio shown by xA quasars, we propose a method to derive the main cosmological parameters based on redshift-independent "virial luminosity" estimates from measurements of emission line widths, roughly equivalent to the luminosity estimates based from line width in early and late type galaxies. A major issue related to the cosmological application of the xA quasar luminosity estimates from line widths is the identification of proper emission lines whose broadening is predominantly virial over a wide range of redshift and luminosity. We report on preliminary developments using the AlIIIλ1860 intermediate ionization line and the Hydrogen Balmer line Hβ as virial broadening estimators, and we briefly discuss the perspective of the method based on xA quasars.

INTRODUCTION
Quasars show a rich spectrum of emission lines coming from a side range of ionic species (e.g., Netzer, 2013). Differences in line profiles, line shifts, line intensities hint at diversity in dynamical conditions and ionization levels of the broad line region (BLR; Sulentic et al. 2000). There is an overall consensus that the Hβ and low-ionization lines (LILs) in general are mainly broadened via Doppler effect from emitting gas virial motion, while high-ionization lines (HILs) are more affected by non-virial broadening, probably associated with a disk wind or clumpy outflows (e.g., Marziani et al., 1996;Elvis, 2000;Proga et al., 2000;Richards et al., 2011, for a variety of approaches). The relative prominence of the wind and disk component reflects the different balance between radiation and gravitation forces (Ferland et al., 2009;Marziani et al., 2010;Netzer and Marziani, 2010), and is systematically changing along the quasar main sequence (MS). If the main driver of the stellar H-R diagram can be considered mass, for the quasar main sequence it is the luminosity to black hole mass ratio that apparently matters most (Boroson and Green, 1992;Sulentic et al., 2000;Marziani et al., 2001;Shen and Ho, 2014;Sun and Shen, 2015;Panda et al., 2019). In addition, unlike stars, quasars are not endowed of spherical symmetry but only by axial symmetry: the spin of the black hole and the axis of the accretion disk define axes of symmetry and a preferential plane (i.e., the equatorial plane of the black hole, the inner accretion disk plane). Aspect becomes an important parameters for AGN. The effect is most noticeable in the distinction between type-1 and type-2 (obscured) AGN (e.g., Antonucci, 1993;Urry and Padovani, 1995), but is relevant also for type-1 AGN (unobscured), as the viewing angle (defined as the angle between the disk axis and the line of sight) may range between 0 and 60 degrees.
Quoting Marziani et al. (2001), we can say that . . . the observed correlation . . . can be accounted for if it is primarily driven by the ratio of AGN luminosity to black hole mass (L/M ∝ L/L Edd , Eddington ratio) convolved with source orientation. Several recent reviews papers have dealt with the quasar main sequence (e.g., Sulentic et al., 2008;Sulentic and Marziani, 2015;Marziani et al., 2018b). Here we will summarize briefly the basic concepts and trends (section 2), afterward focusing on a particular class of objects, the strongest FeII emitters along the sequence (section 3). The basis of the virial luminosity estimations are then presented in section 4, also in an aspect-angle dependent form. Preliminary results for cosmology are illustrated in section 5 using Hβ and the AlIIIλ1860 line to achieve a wide redshift coverage.

A MAIN SEQUENCE FOR TYPE-1 (UNOBSCURED) QUASARS
The Main Sequence stems from the definition of the so-called Eigenvector 1 of quasars by a Principal Component Analysis of PG quasars (Boroson and Green, 1992). The MS of quasars is customarily presented in an optical plane defined by the parameter of FeII prominence R FeII (the intensity ratio between the FeII blend at λ4570 and Hβ) and by the full-width half maximum (FWHM) of Hβ (Boroson and Green, 1992;Sulentic et al., 2000), and associated with the anti-correlation between strength of FeIIλ4570 (or [OIII]λ5007 peak intensity) and width of Hβ. FeII maintains the same relative intensity of the emission blends (at least to a first approximation) but FeII intensity with respect to Hβ changes from object to object, and in the strongest FeII emitters, FeII emission can dominate the thermal balance of the BLR (Marinello et al., 2016). Therefore, it does not appear surprising that the FeII dominance parameters is a fundamental parameter in organizing the diversity of type-1 quasars. The FWHM Hβ is related to the velocity field in the low-ionization part of the BLR, and ultimately associated with the orientation angle, black hole mass and Eddington ratio (section 2.2).
Since 1992, the Eigenvector 1 MS has been found in increasingly larger samples, involving more and more extended sets of multifrequency parameters (Sulentic et al., 2007;Zamfir et al., 2010;Shen and Ho, 2014;Wolf et al., 2019). In addition to the parameters defining the optical plane, R FeII and FWHM(Hβ), several multifrequency parameters related to the accretion process and the accompanying outflows are also correlated (see Fraix-Burnet et al., 2017, for an exhaustive list).
The MS allows for the definition of spectral types (Sulentic et al., 2002;Shen and Ho, 2014). In addition Sulentic et al. (2000) introduced a break at 4,000 km s −1 and the distinction between Population A (FWHM Hβ 4,000 km s −1 ) which includes the Narrow Line Seyfert 1 (NLSy1s), and Population B of broader sources. Extreme Population A (xA) sources are the ones toward the high R FeII end of the main sequence (R FeII > 1). The optical spectra of the prototypical Narrow Line Seyfert 1 (NLSy1) I Zw 1 and the Seyfert-1 nucleus NGC 5548 illustrate the changes that occur along the sequence. I Zw 1 (a prototypical xA source) shows a low-ionization appearance, with strong FeII, weak [OIII]λλ4959,5007, and a spiky Hβ broad profile; on the converse the optical spectrum of NGC 5548 suggests a higher degree of ionization, with low FeII emission, strong [OIII]λλ4959,5007. It is expedient to consider the Hydrogen Balmer line Hβ and the CIVλ1549 as prototypical optical and UV LILs and HILs, respectively. Comparing the CIVλ1549 and Hβ lines yields important information: in the spectra of NGC 5548 the two lines show similar profile, while in the ones of I Zw 1 the CIVλ1549 emission is almost fully blue shifted with respect to a predominantly symmetric and unshifted (with respect to rest frame) Hβ profile (e.g., Sulentic et al., 2000;Leighly, 2004;Coatman et al., 2016).
The comparison between LILs and HILs has provided insightful constraints of the BLR at low-z (Marziani et al., 1996), and this is even more true if the comparison of Hβ and CIVλ1549 is carried out at high L (Shen, 2016;Bisogni et al., 2017;Sulentic et al., 2017;Vietri et al., 2018). Perhaps the most remarkable fact is that a LIL-BLR appears to remain basically virialized (Marziani et al., 2009;Sulentic et al., 2017). The CIVλ1549 blueshift depends on luminosity [the median c( 1 2 ) is ≈ 2,600 and 1,800 km s −1 for Pop. A and B radio-quiet (RQ) at high-L against less than 1,000 km s −1 at low-L], but the dependence is not strong, and can be accounted for in the framework of a radiation driven winds. In spite of the extremely high amplitude CIVλ1549 blueshifts of Pop. A quasars at the highest luminosity L 10 47 erg s −1 , the Hβ profile remains (almost) symmetric and unshifted with respect to rest frame. Figure 1 shows the emission line profiles of Hβ, CIVλ1549, and AlIIIλ1860 and SiIII]λ1892 in three quasars of widely different luminosity. LILs in the optical and UV remains symmetric and unshifted. The lines have been decomposed into two main components.  1994; Popović et al., 2002;Kovačević-Dojčinović and Popović, 2015;Adhikari et al., 2016). Following Sulentic et al. (2002), the BC is modeled by a symmetric and unshifted profile (Lorentzian for Pop. A or Gaussian for Pop. B), as it is believed to be associated with a virialized BLR subsystem. The virialized BLR could be defined as the BLR subregion that is in the condition to emit FeII. Given the stratification of ionization in the BLR, the FeII emission is apparently weighted toward outer radii with respect to Hβ , as also supported by the slightly narrower FeII with respect to Hβ (e.g., Sulentic et al., 2004). The regions emitting Hβ BC 1 and FeII are however believed to be largely co-spatial. • The blue shifted component (BLUE). A strong blue excess in Pop. A CIV profiles is obvious, as in some CIV profileslike the one of the xA prototype I Zw1-BLUE dominates the total emission line flux (Marziani et al., 1996;Leighly and Moore, 2004). For BLUE, there is no evidence of a "stable" profile (in the sense of being represented by the same fitting function), and a "stable" profile is not expected, as we are dealing with a component that is likely subject to phenomena associated with turbulence and hydrodynamic instabilities. In our recent work (Sulentic et al., 2017), we have modeled the blueshifted excess adding to the symmetric BC profile one or more skew-normal distributions (Azzalini and Regoli, 2012). The "asymmetric Gaussian" fitting function is, at present, empirically motivated right by the often-irregular BLUE profile in CIV and MgIIλ2800.
The line width of the LILs clearly increases with luminosity passing from ∼ 1, 000 to 5,000 km s −1 from log L ∼ 44 to ∼ 47. This behavior-typical of Pop. A quasars-suggests that a virialized sub-system emitting mainly LILs coexists with outflowing gas. Even at the highest luminosity the data (Coatman et al., 2016;Sulentic et al., 2017;Vietri et al., 2018) remain consistent with a scenario that posits a dichotomy within the BLR, where the blue shifted emission is produced in clumps of gas or a continuous wind (Collin-Souffrin et al., 1988;Elvis, 2000), and LILs are emitted by gas in virial motion. In this scenario, we observe a net blueshift in CIVλ1549 because of the large gas velocity field component along the radial direction, with the receding side of the flow obscured by an optically thick equatorial structure, most likely associated with the accretion disk. I Zw 1 fits very well this interpretation: if LILs are emitted in the flattened disk structure, the Hβ should appear narrow, while CIVλ1549 emission broader and shifted to the blue, which is indeed the case (Marziani et al., 1996;Leighly, 2004).

The Definition of Spectral Types
The data point occupation in the plane R FeII -FWHM(Hβ) make it expedient to define spectral types; A1, A2... in order of increasing R FeII ; A1, B1, B1+.. . . in order of increasing FWHM(Hβ). This has the considerable advantage that a composite over all spectra within each bin should be representative of objects in similar dynamical and physical conditions. In principle, a prototypical object can be defined for each spectral type to analyze systematic changes along the quasar MS. Before summarizing some basic trends, it is helpful to recall the main motivations behind the distinction between Population A and B.

Population A and NLSy1s
Why choose a FWHM limit for Hβ at 4,000 km s −1 ? Many studies even in recent times distinguish between NLSy1s (FWHM Hβ 2, 000 km s −1 ) and rest of type 1 AGNs (e.g., Cracco et al., 2016). The basic reason to extend the limit to 4,000 km s −1 is that several important properties of NLSy1s are not distinguishable from "the rest of Population A" in the range 2, 000 km s −1 FWHM(Hβ) 4, 000 km s −1 . NLSy1s are located at the lowend of the distribution of FWHM(Hβ) in the samples of Shen et al. (2011) andZamfir et al. (2010) without any apparent discontinuity. NLSy1s and the rest of Population A show no apparent discontinuity as far as line profiles are concerned (Cracco et al., 2016). Composite Hβ profiles of spectral types along the MS are consistent with a Lorentzian as do the rest of Population A sources. A recent analysis of composite line profiles in narrow FWHM(Hβ) ranges (δ FWHM = 1,000 km s −1 ) confirms that there is no discontinuity at 2,000 km s −1 ; best fits of profiles remain Lorentzian up to at least 4,000 km s −1 . A change in the Hβ profile shape occurs around FWHM(Hβ) = 4,000 km s −1 at the Population A limit, not at the one of NLSy1s (Sulentic et al., 2002). No significant difference between CIVλ1549 centroids at half-maximum of NLSy1s and rest of Pop. A is detected using the HST/FOS data of the Sulentic et al. (2007) sample (Marziani et al., 2018a). In addition NLSy1s span a broad range of spectral types, the R FeII parameter being between almost 0 and 2 in the quasi-totality of sources, as in the case of the rest of Population A.
A complete tracing of the MS at high L is still missing (we consider high-luminosity sources the quasars with bolometric log L 47 [erg/s]): the Hβ spectral range is accessible only with IR spectrometers, and high-luminosity quasars are exceedingly rare also at z 1. The main effect of a systematic increase in black hole mass M BH can be predicted. If the motion in the LIL-BLR is predominantly virial, we can write the central black hole mass M BH as a virial mass : which follows from the application of the virial theorem in case all gravitational potential is associated with a central mass. δv K is the virial velocity module, r BLR the radius of the BLR, G the gravitational constant. Equation (1) can be of use if we can relate δv K to the observed velocity dispersion, represented here by the FWHM of the line profile: via the structure factor (a.k.a. f S form or virial factor) whose definition is given by: If the BLR radius follows a scaling power-law with luminosity (r ∝ L a , Kaspi et al. 2000;Bentz et al. 2013), then Or, equivalently, If a = 0.5, the FWHM grows with M BH 0.25 , meaning a factor of 10 (!) for log L passing from 44 (relatively low luminosity) to 48 (very luminous) quasars. The trend may not be detectable in low-z flux limited samples, but becomes appreciable if sources over a wide range in L are considered. At high M BH , the MS becomes displaced toward higher FWHM values in the optical plane of the MS; the displacement most likely accounts for the wedge shaped appearance of the MS if large samples of SDSS quasars are considered (Shen and Ho, 2014).
If we consider a physical criterion for the distinction between Pop. A and B a limiting Eddington ratio (log L/L Edd ∼ 0.1-0.2), then the separation at fixed L/L Edd becomes luminosity dependent. In addition FWHM(Hβ) is strongly sensitive to viewing angle, via the dependence on θ of the f S (section 4.2). An important consequence is that a strict FWHM limit has no clear physical meaning: its interpretation depends on sample properties. At low and moderate luminosity (log L 47 [erg s −1 ]), the limit at 4,000 km s −1 selects mainly sources with log L/L Edd 0.1-0.2 i.e., above a threshold that is also expected to be relevant in terms of accretion disk structure. However, the selection is not rigorous, since a minority of low L/L Edd sources observed pole-on (and hence with FWHM narrower because of a pure orientation effect) are expected to be below the threshold FWHM at 4,000 km s −1 . Such sources include core-dominated radio-loud quasars whose viewing angle θ should be relatively small (Marziani et al., 2001;Zamfir et al., 2008). They are expected to be rare, because the probability of observing a source at a viewing angle θ is P(θ ) ∝ sin θ , although their numbers is increased in flux-limited sample because of the Malmquist effect. As a conclusion, we can state that separating Pop. A and B at 4,000 km s −1 makes sense for low z sample and that, by a fortunate circumstance, Pop. A includes the wide majority of relatively high L/L Edd sources.

EXTREME POPULATION A: SELECTION AND PHYSICAL CONDITIONS
As mentioned in section 1, Eddington ratio is a parameter that is probably behind the R FeII sequence in the optical plane of the MS. It is still unclear why it is so. At a first glance, the low-ionization appearance of sources whose Eddington ratio is expected to be higher is puzzling. Marziani et al. (2001) accounted for this results including several observational trends that could lower the ionization parameter to increase R FeII . More recent results include a careful assessment of the role of metallicity (Panda et al., 2018(Panda et al., , 2019. There is empirical evidence of a correlation between R FeII and L/L Edd , but the issue is compounded with the strong effect of orientation on the line broadening, and hence on M BH and L/L Edd computation. As an estimate of the viewing angle is not possible for individual radio-quiet sources (even if some recent works provide new perspectives), conventional M BH (Equation 1) and L/L Edd estimates following the definitions in the previous section suffer of a large amplitude scatter (Marziani et al., 2019, and references therein). Mass estimates over a broad range of FWHM are especially uncertain, and may lead to sample-dependent, spurious trends. Even if a combination of effects related to Eddington ratio, viewing angle, and metal content reproduces the sequence, and accounts for the relative distribution of sources in each spectral bin, the actual structure of the BLR and its relation with the accretion mode are still unclear. This said, if our basic interpretation of R FeII as mainly related to L/L Edd is correct, then the strongest FeII emitters should be the highest radiators. For instance, Sun and Shen (2015) provide evidence not dependent on the viewing angle: in narrow bins of luminosity, the stellar velocity dispersion of the host spheroid (a proxy for M BH ) is anti correlated with R FeII , implying that L/L Edd increases with FeII. The fundamental plane of accreting black holes (Du et al., 2016) correlates L/L Edd and dimensionless accretion rate with R FeII and with the "Gaussianity" parameter D = FWHM/σ of Hβ. Both methods imply that the highest R FeII corresponds to the highest L/L Edd .
If we select spectral types along the sequence satisfying the condition R FeII 1 (A3, A4, . . . ), we select spectra with distinctively strong FeII emission and Lorentzian Balmer line profiles. They account for ∼ 10% of quasars in low-z, optically selected sample FeII in Pop. A.
The simple selection criterion • R FeII = FeIIλ4570 blend/Hβ > 1.0 corresponds to the selection criterion: The R FeII >1 and the UV selection criteria are believed to be equivalent. Marziani and Sulentic (2014a) compared the distribution of sources satisfying R FeII 1 in the plane AlIII λ1860/SiIII]λ1892 vs. SiIII]λ1892/ CIII]λ1909, and found that these sources were confined in a box satisfying the boundary conditions AlIII λ1860/SiIII]λ1892>0.5 & SiIII]λ1892/ CIII]λ1909>1, as expected if the two selection criteria are consistent. The number of sources for which this check could be carried out is however small, and several xA sources are located in borderline positions. Work is in progress to test the consistency on a larger sample.
Though, the UV xA spectrum is very easily recognizable, as shown by the composite spectrum of Martínez-Aldama et al. (2018). Lines have low equivalent width: some xAs are weak lined quasars (W(CIVλ1549) ≤ 10 Å, WLQ Diamond-Stanic et al., 2009), whereas WLQs can be considered as the extreme of extreme Pop. A. The CIIIλ1909 emission almost disappears. In the plane log U−log n H defined by CLOUDY simulation, UV line intensity ratio converge toward extreme values for density (high, n H > 10 12 −10 13 cm 3 ), and ionization (low, ionization parameter U ∼ 10 −3 − 10 −2.5 ). Extreme values of metallicity are also derived from the intensity ratios CIV/AlIII CIV/HeII AlIII/SiIII] (Negrete et al., 2012;Martínez-Aldama et al., 2018;Sniegowska et al., in preparation), most likely around 20 times solar or with abundances anomalies that increase selectively aluminum or silicon, or both.
In addition to the selection criteria, a virial broadening estimator equivalent to Hβ should be defined from the emission line in the rest frame UV spectrum. The CIVλ1549 line is unsuitable unless heavy corrections are applied (Coatman et al., 2017;Marziani et al., 2019). However, the UV spectrum of xA quasars at z ∼ 2 show symmetric low-ionization and blueshifed high-ionization lines even at the highest luminosity The Eddington ratio precise values depend on the normalization applied; the relevant result is that xA quasars radiate at extreme L/L Edd along the MS with small dispersion (Marziani and Sulentic, 2014a). This results is consistent from the expectation of accretion disks at very large accretion rates. Accretion disk theory predicts low radiative efficiency at high accretion rate and that L/L Edd saturates toward a limiting values (Abramowicz et al., 1988;Mineshige et al., 2000;Sadowski et al., 2014). The L/L Edd distribution shown by Marziani and Sulentic (2014a) provides some support to this finding. While R FeII and L/L Edd might be correlated, as mentioned above, if R FeII 1, L/L Edd apparently scatters around a well-defined value with a relatively small scatter ≈ 0.13 dex. Clearly, this results should be tested by larger sample and, even more importantly, by Eddington ratio estimators not employing the FWHM of the line used for the M BH computation. Another important fact is the self similarity of the spectra selected by the R FeII criterion: as Figure 1 shows, the LILs broadens from 1,000 to over 5,000 km/s, but the relative intensity ratios (and so the overly appearance of the spectrum) remains basically unchanged. This occur over an extremely wide luminosity range, log L ∼ 44 − 48 [erg s −1 ], covering local type 1 quasars in the local Universe as well as the most luminous quasars that are now extinct but that were shining bright at redshift 2 when the Universe had one quarter of its present age.
Accretion disk theory predicts that at high accretion rate a geometrically thick, advection dominated disk should develop (Abramowicz et al., 1988;Sadowski et al., 2014), especially at the radii closest to the central black hole where the temperature should be so high that radiation pressure dominates over gas pressure. The disk vertical structure and the interplay between BLR and disk remains to be clarified (Wang et al., 2014c). It might be one of the biggest challenges for the next decade. The issue can be addressed by two dimensional reverberation mapping and by careful modeling of the coupling between dynamical and physical conditions (Li et al., 2013(Li et al., , 2018Pancoast et al., 2014). A tentative model of xA sources is provided by the sketch of Figure  6 of Marziani et al. (2018b). The innermost part of the disk is puffed up by radiation pressure, while the outermost one remains geometrically thin. This change from the standard thin α disk provides two key elements for the BLR structure: the existence of a collimated cone like region, where the high ionization outflows might be produced. While a vertical shape of the outer side of the geometrically thick region is unlikely to be appropriate, shadowing of the low-ionization emitting region is expected to occur, with the unpleasant consequence that the continuum seen by the BLR may not be the one seen by the observer (Wang et al., 2014c).

EDDINGTON STANDARD CANDLES AS DISTANCE INDICATORS FOR COSMOLOGY
From the discussion of the previous sections we gather that three conditions are satisfied for xA quasars: 1. Constant Eddington ratio L/L Edd ; xA quasars radiate close to Eddington limit: Note that the precise value or L/L Edd depends on the normalization applied for M BH and on the bolometric correction to compute the bolometric luminosity. Applying widely employed scaling laws, L/L Edd 1. 2. The assumption of virial motions of the low-ionization BLR, so that the black hole mass M BH can be expressed by the virial relation (Equation 1): 3. Spectral invariance: for extreme Population A, the ionization parameter U can be written as Netzer (2013), where Q(H) is the number of hydrogen ionizing photons. U has to be approximately constant, otherwise we would observe a significant change in the spectral appearance.

The Virial Luminosity Equation
Taking the three constraints into account, the virial luminosity equation derived by Marziani and Sulentic (2014a) can be written in the form: where the energy value has been normalized to 100 eV (ν i,100eV ≈ 2.42 · 10 16 Hz), κ i,0.5 is the fraction of bolometric luminosity belonging to the ionizing continuum scaled to 0.5, the product (n H U) has been scaled to the typical value 10 9.6 cm −3 (Padovani and Rafanelli, 1988;Matsuoka et al., 2008;Negrete et al., 2012), and the FWHM of the Hβ BC is expressed in units of 1,000 km s −1 . The f S is scaled to the value 2 following the determination of Collin et al. (2006). The FWHM of Hβ broad component and of AlIIIλ1860 are hereafter adopted as a virial broadening estimator. Equation (8) implies that, by a simple measurement of the FWHM of a LIL, we can derive a z−independent estimate of the accretion luminosity Sulentic, 2014a, c.f. Teerikorpi, 2011). Recent works involved similar ideas concerning the possibility of "virial luminosity" estimates, and confirm that extremely accreting quasars could provide suitable distance indicators because their emission properties appear to be stable with their luminosity scaling with black hole mass at a fixed ratio (Wang et al., , 2014aLa Franca et al., 2014). Equation (8) can be applied to xA quasars only (L/L Edd ∼ 1). However, it can be applied to all xA distributed over a wide range of luminosity and redshift, where conventional cosmological distance indicators are not available (see D'Onofrio and Burigana, 2009;Hook, 2013;Czerny et al., 2018, for reviews concerning distance indicators beyond z 1).
The virial luminosity equation is analogous to the Tully-Fisher and the early formulation of the Faber Jackson laws for earlytype galaxies (Faber and Jackson, 1976;Tully and Fisher, 1977). We note in passing that galaxies and even clusters of galaxies are virialized systems that show an overall consistency with a law L ∝ σ n , where σ is a measurement of the velocity dispersion and n ∼ 2 − 4. Both the Faber and Jackson (1976) and Tully and Fisher (1977) found some application to the determination of the distance scale, especially at low-z (e.g., Dressler et al., 1987, see Czerny et al. 2018 for a recent review). Distance measurements to early-type galaxies can be obtained by applying the relation of the fundamental plane of galaxies (e.g., Saulder et al., 2019) which consider the relation between σ , structure and luminosity. Galaxies are more complex than quasars due to their nonhomology . We hope to have minimized structural and physical condition differences by restricting the selection to xA quasars. However, the factor L 0 is most likely affected by intrinsic variations in the spectral energy distribution (SED) of xA sources, as well as by anisotropic emission expected from a thick accretion disk (Wang et al., 2014c). The issue will be briefly discussed in section 6.

Orientation-Dependent Virial Luminosity Equation
The effect of orientation can be quantified by assuming that the line broadening is due to an isotropic component + a flattened component whose velocity field projection along the line of sight is ∝ sin θ : If we considered a flattened distribution of clouds with an isotropic δv iso and a velocity component associated with a rotating flat disk δv K , the structure factor appearing in Equation (8) can be written as which can reach values 1 if κ = δv iso δv K ≪ 1, and if θ is also small ( 30 deg). The assumption f S = 2 implies that we are seeing a highly flattened system (if all parameters in Equation (8) are set to their appropriate values); an isotropic velocity field would yield f S = 0.75.
The virial luminosity equation may be rewritten in the form: where L vir is the true virial luminosity (which implies f S = 1) with L • 0 = 1 4 L 0 , since L 0 was scaled to f S = 2.

The Basic
Following Weinberg (1972Weinberg ( , 2008 the luminosity distance d L can be written as: where M is the energy density associated to matter and the energy density associated to , and d H ≡ c H 0 = 3000 h −1 Mpc = 9.26 × 10 25 h −1 m, with log d H ≈ 28.12[cm], for H 0 = 70 km s −1 Mpc −1 . If M + > 1, S(x) is defined as sin(x) and κ = 1 − M − ; for M + < 1, S(x) = sinh(x) and κ as above; and for M + = 1, S(x) = x and κ = 1. The equation of Perlmutter et al. (1997) has been obtained by posing k + + M = 1, where k is the energy density associated with the curvature of space-time. The luminosity distance d L = d C · (1 + z) can be computed in practice from: with log d H ≈ 28.12 [cm] and The distance modulus µ is defined by that is µ = log d L (H 0 , M , k , , z) + 5, where the constant +5 becomes +25 if distances are expressed in Mpc (Lang, 1980). The distance modulus µ computed from the virial equation yielding L(δv) can be written as: where the constant −2.5 log(4πδ 2 10pc ) =-100.19, with δ 10pc ≈ 3.08 · 10 19 the distance of 10pc expressed in cm. The f λ λ can be the flux at 5100 Å for the Hβ sample, or the flux at 1700 Å if the δv=FWHM comes from the AlIIIλ1860 and SiIII]λ1892 lines. The bolometric correction BC has also to be changed accordingly.
Note that (f λ λ) refers to the rest frame fluxes (the term (1 + z) 2 appears in both sides of the equations, for the virial luminosity and for the distance modulus where the luminosity distance was used).

Results
The Hubble diagram of Figure 3 includes sources from Marziani and Sulentic (2014a), with both Hβ and AlIIIλ1860, and is supplemented by SDSS low-luminosity xA sources at low z covering the Hβ range , by new highz quasar observations at redshift 2 z < 3 that were observed with the OSIRIS camera at the GTC, covering the λ1900 Å blend and hence AlIIIλ1860 (Martínez-Aldama et al., 2018). An additional sample at very low-z is the one of the xA sources considered for the Du et al. (2016) fundamental plane. Preliminary results of the AlIIIλ1860 measurements for ≈ 15 quasars (Sniegowska et al., in preparation) are also shown. The value of log L 0 ≈ 45.06 was used in the computation of the distance modulus µ for this work. The plots in Figure 3 involves a total of 253 sources and indicate a scatter δµ ≈ 1.2 mag. A similar result was shown with 169 xA quasars by Marziani et al. (2017). The plot of Figure 3 shows that even our heterogeneous sample joining Hβ and AlIIIλ1860 based sub-samples already rules out extreme Universes such as the Einstein-de Sitter Universe ( M = 1, = 0, red line in Figure 3), which under-predicts µ values at high z by ≈ 1 mag with respect to concordance. The inclusion of the AlIIIλ1860 has been instrumental to the extension of the Hubble diagram to redshifts 1. Observations of Hβ at redshift 1 are available for several hundreds of quasars (e.g., Coatman et al., 2016;Shen, 2016); this means that only a few tens of spectra may be available from literature, if the observations are restricted to xA source and high S/N is requested. An additional virial broadening estimator is needed. The AlIIIλ1860 line can be measured from SDSS survey up to z 4. Figure 2 shows preliminary results for a sample of Population A sources. The evidence collected so far suggests that the correlation should remain valid for most xA sources (a few notable exceptions are known; for example HE0359-3959, Martínez-Aldama et al. 2017). The suitability of AlIIIλ1860 as a virial broadening estimator will be further explored in eventual works.
The Hubble diagram for quasars restricted to the 92 sources of Marziani and Sulentic (2014a) was also consistent with concordance CDM provided some constraints on M : 0.19 +0.16 −0.08 . The redshift range 2-3 is highly sensitive to M . The new sample is consistent with the old one and allows to set more stringent limits to M ≈ 0.30 +0.06 −0.06 via a standard Bayesian least square fit, employing a sample with no restriction on data (i.e., including low-S/N AlIIIλ1860 spectra and a few quasars at z in the range 2.6-3 not shown in Figure 3).

DISCUSSION
In this paper we have described the method developed by Marziani and Sulentic (2014a) in more detail, presented an expanded Hubble diagram with more objects with respect to the one presented by Marziani et al. (2017), as well as an updated (but still preliminary estimate) of M and the Hubble diagram with concordance and Einstein-de Sitter cosmologies, just to illustrate the sensitivity of the method to cosmological parameter changes. It should be also noted that a comparison between the constraints set by the supernova photometric survey described by Campbell et al. (2013) and a mock sample of 400 quasars with rms = 0.3 uniformly covering the z range 0.1 − 3 that assumes concordance cosmology (shaded contours in Figure 4 of Marziani and Sulentic 2014b) emphasizes the potential ability of the quasar sample to better constrain M .
An elementary error budget suggests that the main uncertainty is associated with FWHM measurement errors and with uncertainty in the structure factor (Marziani and Sulentic, 2014a;Negrete et al., 2018). The broad profile of both LILs and HILs in each quasar spectrum can be modeled by considering two main components, the Hβ BC and the blueshifted component BLUE (section 2). Results shown in Figure 3 refer to Hβ BC profiles whose FWHM has been corrected for the effect of BLUE via a multi-component fit, as described in Negrete et al. (2018). The Hβ BLUE affects however the Hβ FWHM only at ≈ 10% level, with an effect on L that should be 0.3 dex. Vetting the sample on strongest FeII and weakest [OIII]λλ4959,5007 emitters reduces the scatter but only slightly . Negrete et al. (2018) have convincingly shown that orientation effects are at the origin of a large fraction of the scatter if lines are emitted in a flattened system (κ 0.2). The difference between virial and concordance luminosity estimates can be zeroed if the structure factor is assumed to be dependent on angle θ in the form given by Equation (10). The resulting distribution of θ values ranges from 0 to 50 with a peak around θ ≈ 18. The range of viewing angles in the sample of Negrete et al. (2018) is actually modest. The effect of anisotropy from a geometrically and optically thick disk on the SED is not extremely strong within θ 30, approximately around 0.2 dex, following the models of anisotropic emission of Wang et al. (2014c). The actual effect on virial luminosity estimates could be even less than that, as the dispersion in the distribution of θ in the sample of Negrete et al. (2018) is just 7 degrees, which implies that most sources are seen in a viewing angle range between 10 and 30 degrees. Figure 4 of Wang et al. (2014c) indicates that different θ should lead to an increase of the dispersion of a few 0.01 dex around the luminosity expected for 20 degrees at optical wavelengths. This will mainly affect the continuum λf λ that we receive entering in Equation (17).
The presence of a thick disk implies self-shadowing effects, and that the line emitting gas might not be exposed to the same ionizing continuum that is seen by the observer. This second effect related to the disk anisotropy is especially severe if the line emitting gas is located close to the equatorial plane of the disk: in this case, following Wang et al. (2014c), we expect a drastic reduction of the number of ionizing photons which would affect L 0 . Precisely quantifying the effect on line emission is not trivial, as it depends on the distribution of the line emitting gas around the equatorial plane of the disk, in terms of both angular and radial distance (Du and Wang, 2019), the black hole mass as well as the accretion disk model. Differences in the accretion disk SED are also expected from object to object because of different black hole spins (Wang et al., 2014b). The spin and the anisotropy effects on L 0 should be considered in an eventual work, but there is a basic point that we want to remark about the usability of Equation (8) for cosmology.
The SED parameters entering in the relation for the virial luminosity areν i , κ i , and more indirectly, n H U. The ionizing photon flux n H U shows remarkably small changes along the MS (i.e., even without restricting the attention to xA sources), as shown by the various estimates reported in Table 4 of Negrete et al. (2012): it cannot change by a large factor, otherwise the emission line spectrum would change dramatically (Panda et al., 2018(Panda et al., , 2019. The other SED-related parameters appear as κ i /ν i , and are qualitatively affected in the same way by a change of the ionizing continuum. Explorative computations using CLOUDY (Ferland et al., 2017) indicate that the ratio κ i /ν i decreases by a factor ≈ 1.2 for a change of θ from θ 20 − 30 to θ ∼ 50 − 70. This condition could be appropriate if the "walls" of the thick disk are illuminating the low-ionization part of the BLR. Object-byobject differences may add scatter in L 0 , therefore contribution to the dispersion observed in the Hubble diagram of Figure 3.
The issue of anisotropy and of other factors affecting the L 0 term multiplying the line width in Equation (8) brings to the attention the possibility of structural changes with redshift and luminosity although it is at least conceivable that the accretion mode at high rates is giving rise to a stable and reproducible structure. This consideration is supported by the stability of xA sources at least in the optical domain , as shown by the small amplitude of continuum variations detected in dedicated reverberation mapping campaigns, and by the similarity of low-ionization emission line profiles over a wide range of redshift and luminosity.
Summing up, we can say that the proper calibration of the data for cosmology will need large samples of quasar spectra covering the rest frame optical and UV, at high S/N, over a wide range of redshift. This is a conditio sine qua non to test for systematic sources of errors and hence for the effective exploitation of quasar data for cosmology applying the virial luminosity equation at z 1.

CONCLUSION
The quasar MS has allowed us to isolate sources which are the highest radiators among quasars. This is already an important feat, as xAs can then be selected by the application of simple criteria based on optical and UV emission line ratios. xA quasars show a relatively high prevalence (10%) and are easily recognizable in the redshift range 0-5. The present work illustrated several key aspects of the method, and provided the virial luminosity equation dependent from the viewing angle in an explicit form (a result of Negrete et al. 2018). However, several additional caveats need to addressed in subsequent work: the reliability of AlIIIλ1860 as a virial broadening estimator and, generally speaking, the inter-calibration of the data obtained in the rest frame UV and optical domain, as well as the effect of anisotropic continuum emission expected for xA sources. xA might be suitable as Eddington standard candles especially if orientation and systematic effects as a function of redshift and luminosity can be identified and accounted for.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.