The phylogeny of quasars and the ontogeny of their central black holes

The connection between multifrequency quasar observational and physical parameters related to accretion processes is still open to debate. In the last 20 year, Eigenvector 1-based approaches developed since the early papers by Boroson and Green (1992) and Sulentic et al. (2000b) have been proven to be a remarkably powerful tool to investigate this issue, and have led to the definition of a quasar"main sequence". In this paper we perform a cladistic analysis on two samples of 215 and 85 low-z quasars (z 0.7) which were studied in several previous works and which offer a satisfactory coverage of the Eigenvector 1-derived main sequence. The data encompass accurate measurements of observational parameters which represent key aspects associated with the structural diversity of quasars. Cladistics is able to group sources radiating at higher Eddington ratios, as well as to separate radio-quiet (RQ) and radio-loud (RL) quasars. The analysis suggests a black hole mass threshold for powerful radio emission and also properly distinguishes core-dominated and lobe-dominated quasars, in accordance with the basic tenet of RL unification schemes. Considering that black hole mass provides a sort of"arrow of time"of nuclear activity, a phylogenetic interpretation becomes possible if cladistic trees are rooted on black hole mass: the ontogeny of black holes is represented by their monotonic increase in mass. More massive radio-quiet Population B sources at low-z become a more evolved counterpart of Population A i.e., wind dominated sources to which the"local"Narrow-Line Seyfert 1s belong.


INTRODUCTION
Quasars represent the most luminous stable sources in the Universe (e.g., D'Onofrio et al., 2012D'Onofrio et al., , 2016, and it is disturbing that some basic issues still remain unsolved in spite of decade-long efforts (e.g., Sulentic et al., 2012;Antonucci, 2013, and references therein). The overall most pressing issue is perhaps the inability to connect physical parameters (such as black hole mass; Marziani and Sulentic 2012;Shen 2013;Peterson 2014, for recent reviews) to the observational properties of individual quasars. This involves a comprehensive structural and dynamical modeling of the broad line emitting regions (e.g., Netzer, 2013, Fig. 1 shows the main sequence (MS) in the plane Fe II opt prominence (R FeII = I(FeIIλ4570)/I(Hβ)) and FWHM of broad Hβ Sulentic et al., 2000a;Shen and Ho, 2014). Several radio, IR, optical, UV, and X properties of low-z quasars (listed in Table 1) can be organized along the MS. To focus on the main physical aspects, Sulentic et al. (2000b) defined a 4D parameter spaces with parameters that are observationally orthogonal: in addition to R FeII and FWHM(Hβ) which are indicators of the low-ionization emitting gas physical conditions and dynamics, they considered the high-ionization line CIVλ1549 shift with respect to the systemic reference of quasars (a tracer of powerful outflows, Richards et al. 2011), and the soft X-ray photon index Γ S , related to the accretion mode (e.g., Pounds et al., 1995;Grupe, 2004).
The data in Table 1 are reported for typical Population A (FWHM 4000 km s −1 ) and Population B (broader) sources following the definition of Sulentic et al. (2000a). The separation into two populations was originally proposed on the basis of a rather abrupt discontinuity in the shape of the Balmer line profiles (e.g., Sulentic et al., 2002;Collin et al., 2006). Systematic changes at FWHM ∼ 4000 km s −1 are probably associated with a discontinuity in accretion mode , as the most relevant physical parameter governing the MS and the differences summarized in Table 1 is likely Eddington ratio Shen and Ho, 2014).
Observational differences for sources belonging to the two Populations become obvious looking at the extremes of the MS. At the top left positions on the sequence shown in Fig. 1 (Population B), lines are very broad (often with composite and couple peaked profiles, Strateva et al. 2003), redward asymmetric, Fe II opt is weak, the spectral energy distribution (SED) is hard (Laor et al., 1997;Shang et al., 2007), [OIII]λλ4959,5007 is strong and symmetric with a blueward asymmetry close to the line base. At the other end, extreme Pop. A sources (R FeII 1) are high accretors which show evidence of strong radiation driven winds Richards et al., 2011;Coatman et al., 2016).
All these developments in the study of quasars strongly indicate that statistical multivariate analyses could be useful. They should bring a more powerful capacity to determine the relevant properties that explain the observed diversity. In this paper, since quasars are known to be evolving objects, we have chosen to use a phylogenetic approach which establishes evolutionary relationships. Among the possible tools, cladistics is certainly the most general and simplest to implement.
In the following, we first describe the cladistic method ( § 2). The method is then applied to a set of optical and UV parameters that include the four parameters of the Sulentic et al. (2000b) 4D space, as well as several additional key parameters ( § 2.3). They are available for a sample of 85 sources (Sect. 2.2) which is as subsample of 215 sources for which only the optical spectral range is available and that is also considered in this study. The cladistic analysis resolves the difference between Pop. A and B and indicates that RL sources are the most evolved Population B sources ( § 3). The results are briefly discusses in terms of quasar populations at low-z as well as of evolution over cosmic age ( §4.2), and in light of our present understanding of the RQ/RL dichotomy ( §4.3).

Cladistic Analysis
Astrocladistics aims at introducing phylogenetic tools in astrophysics. A phylogeny shows the evolutionary relationships between groups, species or classes, while a genealogy shows the relationships between individuals with parents and offsprings. Ontogeny is the evolution of a single individual. In astrophysics, the observations only allow for phylogenetic reconstructions. Consequently, in this paper, each quasar supposedly represents a species (or sub-species, class...). Nevertheless, for simple objects, like a black hole or a star, it is possible to observe identical objects at different stages of evolution, so that their ontogeny can be derived. Among the tools developed for phylogenetic analyses, cladistics, also called Maximum Parsimony, is the most general and the simplest to implement. It uses parameters, and not distances, to establish the relationships between the species by minimizing the total evolutionary cost depicted on a phylogenetic tree.

Outline of the astrocladistics approach
The cladistic analysis is a phylogenetic method for classification, an unsupervised multivariate classification technique that establishes the relationships between the taxa under study. A taxon is a class, a species, or an individual supposedly representing a class. The relationships are depicted on a phylogenetic tree that represents the simplest evolutionary scenario given the data.
Contrarily to many clustering and phylogenetic techniques, cladistics does not compute distances between the taxa, but uses the parameters themselves. These parameters ideally must bear an information of evolution, so that they can be discretized into evolutionary stages. In this case they are called "characters", and the number of changes in stages (steps) represents an evolutionary cost. For a given tree, the total number of steps, considering all the characters and all the taxa, defines the complexity of the evolutionary scenario depicted by the tree. The cladistic algorithm aims at finding the simplest evolutionary scenario which is given by the most parsimonious tree, the one that has the lowest total number of steps, among all possible tree topologies that can be constructed with the taxa in the sample. This kind of tree is called a cladogram. Because it is not based on distances, the cladistic algorithm accepts undocumented values.
The reader is referred to Fraix-Burnet et al. (2015) for a review on unsupervised classification techniques in extragalactic astronomy, and to Fraix-Burnet et al. (2006a,b, 2012; Fraix-Burnet (2016) for more details on the use of cladistics in astrophysics. More explanations are available at https://astrocladistics.org.

Implementation of the cladistic analyses
As usual in astrocladistics, we discretized each parameter in each sample into 32 bins in order to keep its continuous and quantative nature. The optimization criterion that we use considers that the evolutionary cost between two states is the absolute value of their difference (l1-norm). The computation used the heuristic search algorithm implemented in Phylogenetic Analysis Using Parsimony PAUP*4.0b10 (Swofford, 2003), with a ratchet method to avoid as much as possible local minima (Nixon, 1999).

Reliability assessment
The reliability ("robustness") of the cladistic analysis has been estimated through several complementary analyses taking slightly different subsets of quasars and parameters. The two samples presented in this paper are the best ones for two reasons. Firstly, the consensus tree obtained from all the most parsimonious trees found for each sample is nearly entirely resolved, showing that they are all in excellent agreement. This indicates that the tree structure corresponds very probably to a global minimum. Secondly, the results with these two samples agree very well with the complementary analyses, some slight disagreement occurring at the relative placement of some groups on the trees. This does not affect the main interpretation of the proposed phylogeny of quasars.

Interpretation of the tree
The phylogenetic tree depicts the shortest path to relate all the objects of the sample through an evolutionary history. Starting from any taxon (a leave on the tree), one can estimate the cost to transform this object into any other one.
The tree is a hierarchical organization of the taxa, and there is no objective way to define species or classes. The substructures of the tree are a good indication that the related objects may be close from an evolutionary point of view. Ideally, a clade is composed of all the taxa that are the descendants emerging from an internal node (i.e. an unlabelled node). These taxa are supposed to have inherited a property from a common (unknown) ancestor placed at the internal node. Hence each substructure has some chance to correspond to a clade and are used here to define the evolutionary groups of quasars.
The detailed structure of the tree depends somewhat on whether the sense of evolution has been imposed. This is made by defining a root, i.e. a taxon or a clade, which is the closest to the ancestor common to all the object of the sample. This is of course not an easy choice to make in a multivariate pattern. One may use a few parameters that are known to evolve in a monotonic fashion, such as the mass of the black hole which cannot decrease easily if at all. One must be careful that this approach does not lead to a conflict with other parameters.
In any case, once the tree is built and can be trusted statistically, it should be seen as a phylogenetic hypothesis extracted from the data, an hypothesis to be confronted with the current astrophysical knowledge.

Sample selection
The selected samples have the non-negligible advantages to cover relatively well the spectral bins of largest occupation along the E1 MS ( Fig. 1), i.e., to provide a fair representation of the quasar spectroscopic diversity at low-z. The evolution of the quasar luminosity function derived by Boyle et al. (2000) for z 0.7 is modest. Quasars with R FeII 1.5 are 98 % of all low-z quasars (Marziani et al., 2013a). It is worth noting that both samples are strongly biased in favor of RL quasars. This is not really a hindrance since an optically-selected, flux-limited complete sample of 85 sources should include only 6 -7 RL sources, encompassing both core-dominated (CD) and lobe-dominated (LD) quasars. The main results related to the MS have been confirmed by the eventual analysis of large Sloan Digital Sky Survey (SDSS)-based samples (e.g., Zamfir et al., 2010;Shen and Ho, 2014).
The larger sample (M215) includes 215 low-z quasars (z 0.7) presented by Marziani et al. (2003a). Measurements of [OIII]λ5007 are available for most of these sources in addition to FWHM(Hβ) and R FeII but the sample lacks UV and soft-X ray information. A second sample (M85) includes 85 sources, and is obtained from the intersection of the Marziani et al. (2003a, 215 sources covering the Hβ spectral range) and the Sulentic et al. (2007, 130 sources with CIVλ1549 covered from Hubble Space Telescope Faint Object Spectrograph (HST/FOS) observations) samples. The measurements used in this paper are available on Vizier for both Marziani et al. (2003a) 1 and Sulentic et al. (2007) 2 . The two samples provide reliable measurements on high S/N, moderate dispersion spectra that are unmatched by the automated measurements obtained on the SDSS spectra by other authors. M85 is necessary for characterizing the dynamics of the line emitting regions, as there is strong evidence that a virialized component (traced by Hβ) is coexisting with an outflowing component traced by CIVλ1549(e.g., Marziani et al., 2016a, Sulentic et al. 2017. The black hole mass (M BH ) has been estimated from Hβ FWHM using the Vestergaard and Peterson (2006) scaling law, and the Eddington ratio computed assuming a factor 10 bolometric correction for the measured specific flux at 5100Å i.e., applying a correction consistent with past and more recent estimates (Elvis et al., 1994;Richards et al., 2006). Figure 1. The optical plane of the 4DE1, FWHM(Hβ) vs. R FeII for the M215 sample. Sources belonging to M85 are in blue, and radio loud sources are represented by circled symbols. The limit of Pop. A and B are marked, as well as the conventional limit of narrow-line Seyfert 1s (NLSy1).

Parameter selection for the cladistic analysis
In a multivariate analysis such as classification, great care should be taken to avoid redundancies or uninformative parameters that could perturb the result. As already said, the characters are expected to trace out the evolution of the quasars, and obviously, correlated parameters may put too much weight on the underlying physical process. However some correlations are not causal (Fraix-Burnet, 2011) and thus should not be eliminated. The disturbing parameters bring noise and can prevent the convergence of the analysis. It is thus important to understand well the parameters for the sample under study and then repeat the cladistics analysis with several subsets of parameters to test the robustness of the phylogenetic signal we are looking for. For the interpretation of the tree, naturally, all available information can be used and any parameters can be projected onto the tree.
Among the available parameters here, the absolute B magnitude M B and L bol are both measures of luminosity and consequently are strongly correlated. We keep only L bol for the analysis. The black hole mass M BH is a derived property and it is included in the results presented in this paper, but analyses performed without this parameter yields very similar results with a slightly less resolved tree. The Eddington ratio is also a derived quantity and was not included in the cladistic analysis presented in this paper. If we include the Eddington ratio in addition to M BH , then the trees become more linear, pointing to a correlation between two or more parameters. As a consequence, we prefer not to use it to establish the tree, and include it for the interpretation only.
The M85 sample includes eleven parameters, that is four additional ones with respect to M215: the photon index Γ S (Gamma), the radial velocity centroid displacement of CIVλ1549 at half maximum c( 1 2 ) (c1o2CIV), and the rest-frame equivalent width of CIVλ1549, W(CIVλ1549) (WCIV). There are twelve unknown values for voIII and four for Gamma.

RESULTS
The result of a cladistic analysis is an unrooted tree. It is however convenient to root it for the definition of groups and for a possible evolutionary interpretation. In the case of quasars, a parameter that can root the cladistic trees is black hole mass, since M BH can only grow as a function of cosmic time: the only way a black hole can disappear is through emission of Hawking's radiation (Hawking, 1974) which is tremendously inefficient for massive black holes. A tree rooted in this way shows a clustering consistent with evolution from less massive to more massive sources. This evolution would represent the real diversification arrow of quasars if and only if the black hole mass is a reliable evolutionary clock. This would in particular imply that the common ancestor of all quasars has a small black hole. We do not make this assumption in this paper since this parameter is merely one part of the numerous and complex transformation processes of these objects. With this caution stated, the quasar sample contains a population of massive quasars which can be seen as "more evolved" than a population of less-massive quasars that are radiating at a higher L/L Edd . The groups can be defined following the substructures of the tree. Formally, clades are (monophyletic) groups of taxa including an hypothetical ancestor (an internal node) and all its descendants. There is no constraint on the level of granularity of the classification, except for homogeneity of the groups and convenience for the interpretation.
We have chosen the groups as shown in Figs. 2 and 3. We find 20 groups on the M85 tree and 41 on the M215 one. 3 The relative homogeneity of the groups can be estimated from the barplots drawn in front of the trees, as well as some evolutionary trends along this phylogeny. In addition, there are individual terminal branches (or leaves, i.e. leading to only one quasar) that are identified with black color in this paper. Each can be considered as representative of a group, increasing the total number of identified groups. Note that the color progressions are identical in the two cases. The agreement between the two classifications can be estimated by comparing them with the M85 sample objects only, with the Adjusted Rand Index that measures the number of pairs of objects that belong or not to the same class in the two classifications. In a perfect match it should be 1, and 0 for random classifications. Here we find 0.42, that is quite good since the two classifications do not use the same set of parameters. This can also be verified on the boxplots. , W(CIVλ1549), in addition to the seven parameters used to find the tree. Colors correspond to the groups defined from the tree on the left. They do not match in any way the groups and colors for the M85 sample, except for the progression from red to deep blue Each group can be characterized by the statistics of the parameters and shown as boxplots ( Fig. 4 and 5). From this, we can already identify known populations (Fig. 6): in M85, group 1 is due to low-luminosity Pop. A sources, group 2 is dominated by extreme Pop. A, while groups 19 and 20 include almost exclusively radio-loud objects: mainly LDs (20)   Largest R K are associated with the largest M BH , and LD and CD are separated into two distinct groups (here group # 19 and 20 for M85).
At the same time M BH increases quite monotonically, this regularity being independent of the choice to root the tree with the lowest M BH . The trend of M BH along the trees is remarkable: there is a strong increase first and then a kind of plateau that may even be split from the first part. This behavior seems strongly related to the trend in R K , L/L Edd and FWHM(Hβ). The anomalous group 11 in M85 has large R K , W(CIVλ1549): most likely the reflection of large extinction internal to the quasars belonging to this small group. Apart from this, it is apparently not different from the rest of Population B groups.

A threshold M BH for RLness
RL quasars are predominantly found among Pop. B, where RQ quasars are also found with similar accretion parameters (45 log L bol 46 erg s −1 , 8.5 log M BH 9.5, L/L Edd ∼ 0.1). Powerful RL sources appear in our low-z sample only for large M BH , and CD and LD are separated because of an intervening role of orientation (which affects FWHM(Hβ), and L bol ; in phylogenetic terms, CD and LD belong to different monophyletic groups). We insist on the fact that the other analyses that we have performed without M BH give the same results, so that this conclusion does not depend on the presence of this parameter in the analysis nor on its selection for the rooting of the tree. The M BH result is clearly dependent on sample properties. M215 and M85 contains only very powerful radio-loud sources, so that results from additional studies need to be considered to decide whether a mass threshold is appropriate ( §4.3).
The most established tenet of radio-loud unification (Urry and Padovani, 1995) is that the difference between CDs and LDs is due to relativistic beaming and therefore strongly dependent on orientation. It is also known that CD and LD sources show differences in optical spectroscopic properties (e.g.,s Wills and Browne, 1986;Zamfir et al., 2008;Buttiglione et al., 2010). We do not really see a separation that may be ascribed to orientation effects in the other branches and trunk of the cladistic tree. Therefore orientation appears to be an intervening physical parameter that differentiates the RL sources into CDs and LD staring from a common progenitor species.

L/L Edd
Several observational parameters are strongly correlated with L/L Edd Kuraszkiewicz et al., 2004;Shen and Ho, 2014). L/L Edd expresses the ratio between radiation and gravitational forces (Netzer and Marziani, 2010;Marziani et al., 2010). The transition between Pop. A and B occurs at L/L Edd ≈ 0.2 ± 0.1 , which is consistent with the limit that corresponds to the transition from a geometrically thin to a geometrically thick accretion flow sustained by radiation pressure (Abramowicz and Straub, 2014, and references therein). Prominent winds are associated with high L/L Edd Richards et al., 2011;Coatman et al., 2016). Therefore, we can think of two quasars populations, one with relatively modest M BH (7 log M BH 8 [M ]) radiating at high Eddington ratio and associated with strong wind, and one of more massive sources (8 log M BH 10 [M ]) radiating at L/L Edd 0.1. While L/L Edd appears to be the main physical factor governing E1, high-M BH quasars may have resembled low-M BH quasars in an earlier stage of their evolution when they were "wind-dominated" (as further discussed in §4.2).

Population B as evolved Population A quasars
Every quasar is the direct descendant of a seed black hole. However, black holes of masses M BH 10 5 M are extremely difficult to detect if they are placed in the nuclei of external galaxies. Even if they are radiating at L/L Edd ∼ 1 (L/L Edd 2 -3 may not be possible Mineshige et al., 2000), their apparent V magnitude will be ≈ 22 at redshift z ≈ 0.3. In flux-limited quasar samples at low-z, we detect quasars in the mass range 6 log M BH 8 radiating close to their Eddington limit. The masses of these sources are clearly not the masses of the fledgling seed BHs. Nonetheless, in the local Universe, the only sources radiating close to the Eddington limit are these relatively low-M BH quasars. It is easy to see that, if the most massive black holes were nowadays radiating at their maximum radiative power per unit mass (≈ 2 L/L Edd ), they would be almost visible to the naked eye (with M BH ≈ 10 10 M at z ≈ 0.15 it would be m V ≈ 6.7!). The absence of massive BHs radiating close to the Eddington limit is associated with the overall downsizing of the star formation and nuclear activity at recent cosmic epochs (e.g., Fontanot et al., 2009;Reviglio and Helfand, 2009;Hirschmann et al., 2014). The very massive black holes that were shining bright mostly belong now to spent systems (Lynden-Bell, 1969), accreting at a very low rate. In this respect sources like Messier 87 (which we consider a prototypical example of the "spent" quasars which is hosting one of the most massive BHs in the local Universe, Walsh et al. 2013) are very different from the Pop. B sources that we are considering in our sample. Pop. B sources are quasars accreting at a modest pace but high enough to be in a radiatively efficient accretion mode that can be modeled by a geometrically thin, optically thick α-disk with an efficiency η ∼ 0.07 (Shakura and Sunyaev, 1973). Assuming that Pop. B sources are accreting at a constant mass rate, M BH grows in a linear regime, the time needed for a Pop. B source to have grown from a mass M 0 BH = 10 8 M is where f is the ratio of the actual mass to the initial one and we have assumed η 1 (Netzer, 2013, p. 297). The ∆t needed for a BH to grow from the typical masses of the highly accreting quasars at a typical luminosity of luminous low-z quasars is short relative to the expectation of the cosmic evolution of quasar accretion rates. For example, a Pop. A quasar with M 0 BH ≈ 10 8 M at z ≈ 0.2 could be seen as a Pop. B source with M BH ≈ 10 9 M at z ≈ 0.15 since ∆t ≈ 5.1 · 10 8 yr. Its spectroscopic appearance would change as well, moving the source from the bottom right toward the top left of the E1 MS, as schematically indicated in Fig. 7. In other words, it is legitimate to assume that Pop. B are evolved Pop. A sources and not the same quasars that were once radiating at very high L at the cosmic peak of star formation rate and quasar activity.
This eventuality is rather unlikely also if we consider the expected comoving density of spent (dead) quasars in the local Universe. It is convenient to follow the Boyle et al. (2000) luminosity function parameterization: where the evolution is given by the redshift dependence of the break magnitude M * B (z) = M * B (0) − 2.5(k 1 z + k 2 z 2 ), and a ≈ −3.41, b ≈ −1.58. The most luminous quasars are the ones with the most massive black holes radiating close to the Eddington limit (which were relatively rare at z ≈ 2 and have all but disappeared at z 0.7. It is reasonable to assume that all extremely luminous quasars have M BH 10 9.5 M Natarajan and Treister, 2009;King, 2016), and that they may not radiate at highly super-Eddington ratios, the limit being close to a few times the Eddington ratio. For a typical quasar SED, log L ≈ 36.54 − 0.4M B . Therefore, with L Edd ≈ 1.3 · 10 5 M BH M L , M B,max ≈ −27.9 for M BH = 10 9 M and L/L Edd = 1. The comoving number densityñ of sources with M B below this absolute magnitude limit yields a number of sources dN in the comoving volume element dV t, and the number of sources N lum which are supermassive (M BH 10 9 M ) and that were once radiating close to their Eddington limit can be obtained by integrating over 1.5 z 2.3: where d C is the comoving distance, and E(z) = Ω M (1 + z) 3 + Ω Λ . This number can be compared to the number N B of Pop. B sources up to z ≈ 0.75 i.e., to N B = z 2 z 1ñ dV , with z 1 ≈ 0.1 and z 2 ≈ 0.75. Assuming that Pop. B are about 1 2 of all type-1 quasars of all quasars within M B,max ≤ −21 (as suggested by flux-limited optically selected samples Zamfir et al. 2010), the number of very massive black holes that were once radiating above -27.9 in the redshift range 1.5 ≤ z ≤ 2.3 falls short by slightly less than two orders of magnitude to explain the density of the "present-day" (z 0.75) Population B quasars: N lum ≈ 0.03N B . This is hardly surprising since the comoving number density of quasars at M B ≈ −28 and z 1.5 is log Φ(M B ) ≈ −8.15 [Mpc −3 mag −1 ], and at z ≈ 0.5 and M B ≈ −23 is ≈ −6.25 [Mpc −3 mag −1 ]. If the integration limit used to compute N lum from -27.9 is lowered down to M B ≈ -26 at high-z, and if the luminosity function is extrapolated to this absolute magnitude, the number of quasars N lum becomes comparable to N B at low z. Thus, we cannot exclude that a quasar reached M BH ∼ 10 9 M at z ∼ 2, stopped accreting, and was then rejuvenated at recent cosmic times.

From seed to dead supermassive black holes
Recent deep observations of faint quasars (m V ∼ 22) obtained with GTC indicate the presence of a slowly evolving quasar Population at z ≈ 2, not dissimilar to the one observed at low-z as far as the frequency of Pop. B sources is concerned . This result is consistent with the existence of a population of sources that were evolving on timescales much shorter than the Hubble time, then (at a cosmic age of just ≈ 3 Gyr) and now, and with the relatively short lifetimes expected for quasars (Kelly et al., 2010). As mentioned earlier, the comoving density of the very luminous quasars is very low, since they appear at the high end of the luminosity function. Therefore, it may be that -if we exclude the most extreme sources -we are seeing a process going on systematically over cosmic epoch with quasar formation and evolution occurring in a way not so much different to the one expected for a human population (looking back in time, we can still identify adults and young adults as we see them at present). Mathur (2000) and, independently, Sulentic et al. (2000a) suggested that the local-Universe NLSy1 sources accreting at a high rate are reminiscent of the early quasars. We are still far from detecting the first population of "infant" quasars without "adults" (a feat that may never become possible). Unfortunately, at high z it is still not possible to detect black holes of M BH ∼ 10 7 M , even if they are radiating at or slightly above L/L Edd , as there is a redshift-dependent cut-off in the detectable L/L Edd . Our view of quasars is biased. At intermediate to high-z we are "blinded" by the most luminous quasars in the Universe (such as the ones revealed by the Hamburg ESO survey, Wisotzki et al. 2000), objects whose luminosity may also be enhanced by anisotropic emission (Urry et al., 1991;DiPompeo et al., 2014). However, these sources are the ones expected to exert a feedback effect significant enough to affect the host galaxy.
Can we say that the local quasars are repeating the evolutionary patterns of these most luminous quasars? In the context of quasars, feedback effects significant enough to lead to the expulsion of matter on galactic scale may be possible only if the radiative and kinetic powers are extremely high (e.g., Wagner et al., 2013;King and Pounds, 2015, see also the discussions in Marziani et al. 2016a,b). A significant feedback effect of present-day luminous AGNs is proven only in their circumnuclear regions, as shown in the very detailed study of the outflow in NGC 5548 (Kaastra et al., 2014), and energies involved in dispersing the gas outside of the bulge and in a host disk may be beyond the reach of present-day quasars. Why today we do not see them as the more luminous? The answer is related to the lack of available gas to get high accretion rates, and to the dramatic decrease of the expected merger rate with cosmic epoch (e.g., Cavaliere and Vittorini, 2000;Hopkins et al., 2006). In keeping with a human population analogy, we can say that we see a population of young adults and adults which are similar, one close and one at a larger distance. At the larger distance we see a relatively rare population of quasars that are now extinguished. They grew so large that they destroyed their habitat and eventually starved to death.
The cladogram we find shows that if evolution there is, it is very probably along the increasing M BH , with a strong support from the physics and the demography of observed quasars. The difficulty is that the ontogeny describes the evolution of a single quasar. The demonstration in this section shows that a given quasar can only increase M BH and in a relatively short time. If it appears as a Pop. A member, then it will likely evolve toward the Pop. B population. But the cladogram is not supposed to show this; rather, it depicts the phylogeny: if the rooting corresponds to the smaller M BH , then the Pop. A quasars are closer to the most primitive quasars (as far as the diversity in our sample can tell of the diversity of the true population of quasars). This implies several things. First, a Pop. A quasars might not necessarily evolve toward Pop. B, there can be kind of dead-ends if for instance there is no more matter to accrete. These paths could be represented by (some of) the many branches of the trees. Second, and this is the main difference with ontogeny, there are many possible paths from the most ancestral groups. In other words, observing a "young adult" does not tell exactly how the "adult" will look like, even if physically, we can expect that a Pop. A member will have the typical Pop. B appearance. Finally, even if this possibility is probably not imaginable for local quasars, there is no element, from our study, to argue that a quasar cannot appear as such with a very high M BH and high luminosity. This is the limitation of rooting the tree with only one parameter, since quasar appearance cannot be reduced to the mass of their black hole. However, the phylogeny of local quasars as revealed on our trees closely matches the ontogeny of their central black hole. This fact is supported, as mentioned at the end of the previous paragraph, by a physical motivation that emerges from the decrease of gas fraction and merger rate with increasing cosmic epoch.

Radio-loud and radio-quiet quasars
Our cladogram reveals a M BH threshold beyond which radio loudness could be triggered ( § 3.1). This is different from the threshold in the definition of RL quasars for which Zamfir et al. (2008) suggested a value that exceeds by an order of magnitude the canonical limit introduced by Kellermann et al. (1989): log R K 1.8, and/or log P 31.6 erg s −1 Hz −1 . RL sources in our sample are very powerful non thermal emitters and strongly deviate from the correlation between far-infrared (FIR) and radio expected for non-active galaxies (Bonzini et al., 2015;Padovani, 2016). The threshold we find corresponds to the occurrence of RL among the most massive objects, as has been already proposed by several workers (e.g., Chiaberge and Marconi, 2011, and references therein), although this notion is apparently challenged by the discovery of the so-called RL narrow-line Seyfert 1s (e.g., Komossa et al., 2006;Berton et al., 2016).
At the same time, we cannot forget that in the optical plane of E1 we see sources that are both RQ and RL in the same area, in the radiatively-efficient domain associated with a thin accretion disk. Some features such as a strong redward asymmetry in the Hβ profile , are observed in RQ and RL sources alike. RL and RQ Pop. B occupy also the same range in M BH and L/L Edd , with modest L/L Edd , Figure 7. A possible evolutionary scheme for RQ quasars. As M BH increases, sources may move from the location of Population A sources radiating at relatively high Eddington ratio and with much evidence suggesting the presence of a radiation driven wind Sulentic et al., 2007), to Population B sources in the optical plane of the 4DE1 parameter space. as shown in samples larger than the one of the present paper (Sikora et al., 2007). Therefore, assuming that there is a M BH threshold that makes it possible for an AGN to become RL, this condition is necessary but not sufficient. The threshold idea is appealing because it allows for the solution of the MCG 06-30-15 paradox: a maximally rotating black hole inferred from the Fe Kα profile (Iwasawa et al., 1996;Sulentic et al., 1998a) in a radio quiet quasar (actually an anonymous Seyfert 1 with an undistinguished spectrum, Sulentic et al. 1998b). The spin condition is therefore not a sufficient condition to guarantee radio loudness.
In the context of the Blandford-Znajek mechanism for relativistic jet creation (Blandford and Znajek, 1977), the jet power is related to M BH , the spin angular momentum J, and the magnetic field by: where B is the magnetic field and J max = GM 2 /c. This equation indicates that even for a maximally rotating BH, and large mass (M BH ∼ 10 9 M ), a strong magnetic field is needed (∼ 10 2 − 10 4 G).
Recent theoretical work and numerical simulations indicate that jet collimation is occurring in a magnetically-arrested accretion disk (MAD) regime (Punsly, 2015) where magnetic flux is maximized near the black hole, after being dragged from the outer disk toward the center. A modest field strength expected to be present in the nuclei of galaxies (∼ 1µ G) may yield a field strength close to the BH high enough to collimate the jet (Narayan et al., 2003). However, whether a MAD will actually develop, depends on how the magnetic fields associated with the disk behaves as it is dragged toward the central black hole, a still debated issue. It is not understood under which conditions the disk may be able to maintain its magnetization (e.g., Blandford and Payne, 1982;Riols et al., 2016).
Compact steep spectrum and Giga-Hertz peaked sources (O'Dea, 1998) which often show evidence of high L/L Edd from optical and UV properties (Wu, 2009) are accounted for if MAD conditions and jet collimation also occur in an advection-dominated accretion flow (ADAF) context (e.g., Czerny and You, 2016, and references therein) at high accretion rate i.e., for a geometrically thick, optically thick disk. A typical case may be 3C 57, which hosts a very massive BH (M BH ∼ 10 9 M , , is very powerful at present and also shows evidence of relic radio emission, probably due to a past activity cycle (Punsly et al., 2016). Rejuvenated sources may be only a minority of cases, and hence account for the fact the radio-bright population is mostly due to massive black holes radiating at modest L/L Edd , as it is the case of the RL sources in M215 and M85. Therefore, if the attention is restricted to the most powerful RL sources, jet collimation may be made possible by the concomitant occurrence of high mass, high J, strong and ordered magnetic field. A corollary is that the RQ sources in the same L/L Edd and M BH domain should have low values of J or of B. According to recent models, a relatively low L/L Edd ( 0.1) may not be a necessary condition, although it is highly unlikely to have high L/L Edd for very massive BHs at low-z.

CONCLUSION
The cladistic analysis presented a new view of the quasar correlation space based on the first eigenvector of Boroson and Green (1992) and of Sulentic et al. (2000a). In particular, it has been possible: (1) to identify the distinction between Pop. A and B, in a sequence of increasing M BH ; (2) to isolate RL sources among the most massive (more evolved) sources; (3) to separate CD and LD as monophyletic groups having the same progenitor. We infer, from these and previous results, that Pop. B sources may indeed be seen as evolved Pop. A. We tentatively suggest that the radio quiet / radio loud dichotomy is influenced by differences in spin as well as by a different magnetization of the accretion disk.