Didier Fraix-Burnet
email: didier.fraix-burnet@univ-grenoble-alpes.fr

Paola Marziani

Mauro D ' Onofrio

Deborah Dultzin

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

Keywords: black hole physics -quasars, general -quasars, evolution -quasars, radio loud -radio jet formation -multivariate analysis -cladistics

 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., 2012[START_REF] D'onofrio | New Eyes for Galaxies Investigation[END_REF], and it is disturbing that some basic issues still remain unsolved in spite of decade-long efforts (e.g., [START_REF] Sulentic | Fifty Years of Quasars: From Early Observations and Ideas to Future Research[END_REF]Antonucci, 2013, and references therein). The overall most pressing issue is perhaps the inability to connect physical parameters (such as black hole mass; [START_REF] Marziani | Estimating black hole masses in quasars using broad optical and UV emission lines[END_REF][START_REF] Shen | The mass of quasars[END_REF][START_REF] Peterson | Measuring the Masses of Supermassive Black Holes[END_REF], 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., [START_REF] Netzer | The Physics and Evolution of Active Galactic Nuclei[END_REF] and references therein) which is still missing. In addition, there is no clear explanation of the existence of a RL and a RQ quasar population (see e.g., [START_REF] Padovani | The faint radio sky: radio astronomy becomes mainstream[END_REF], for a recent review).

During the last two decades it has become clear that one single model cannot fit all type-1 RQ quasars, and that the diversity of most observational properties cannot be associated only with a viewing angle (i.e., the orientation angle between the accretion disk axis and the line-of-sight). The landmark work of [START_REF] Boroson | The emission-line properties of low-redshift quasi-stellar objects[END_REF] identified two main Eigenvectors in a space of parameters defined by radio, optical, UV and X quasars properties of Palomar-Green quasars. The second eigenvector was found to be correlated with luminosity i.e., a reformulation of the well-known Baldwin effect, the anti-correlation between high-ionization line equivalent width and luminosity [START_REF] Baldwin | Relative quasar luminosities determined from emission line strengths[END_REF][START_REF] Bian | The C IV Baldwin effect in quasistellar objects from Seventh Data Release of the Sloan Digital Sky Survey[END_REF]. The interpretation of the first eigenvector (hereafter E1), dominated by an anti-correlation between full-width at half-maximum FWHM(Hβ) and [OIII]λλ4959,5007 and Fe II opt prominence, was however not so clear since the beginning, and has fostered speculations on the role of black hole mass, chemical composition, black hole spin, Eddington ratio, and orientation (see Sulentic and Marziani, 2015, for a recent review). While orientation effects are probably present in RQ quasars [START_REF] Netzer | Quasar discs. III -Line and continuum correlations[END_REF], they cannot account for the great diversity of observational properties of the line emitting regions, some of them pointing to differences in physical condition, chemical composition and dynamics (e.g., [START_REF] Elvis | A Structure for Quasars[END_REF][START_REF] Marziani | Searching for the Physical Drivers of the Eigenvector 1 Correlation Space[END_REF][START_REF] Chelouche | Continuum shielding and flow dynamics in active galactic nuclei[END_REF][START_REF] Nagao | The evolution of the broad-line region among SDSS quasars[END_REF][START_REF] Netzer | Cosmic Evolution of Mass Accretion Rate and Metallicity in Active Galactic Nuclei[END_REF][START_REF] Shen | The diversity of quasars unified by accretion and orientation[END_REF].

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β [START_REF] Boroson | The emission-line properties of low-redshift quasi-stellar objects[END_REF]Sulentic et al., 2000a;[START_REF] Shen | The diversity of quasars unified by accretion and orientation[END_REF]. 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, [START_REF] Sulentic | The Demise of the Classical Broad-Line Region in the Luminous Quasar PG 1416-129[END_REF] 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, [START_REF] Richards | Unification of Luminous Type 1 Quasars through C IV Emission[END_REF], and the soft X-ray photon index Γ S , related to the accretion mode (e.g., [START_REF] Pounds | RE 1034+39: a high-state Seyfert galaxy?[END_REF][START_REF] Grupe | A Complete Sample of Soft X-Ray-selected AGNs[END_REF].

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., [START_REF] Sulentic | Average Quasar Spectra in the Context of Eigenvector 1[END_REF][START_REF] Collin | Systematic effects in measurement of black hole masses by emission-line reverberation of active galactic nuclei: Eddington ratio and inclination[END_REF]. Systematic changes at FWHM ∼ 4000 km s -1 are probably associated with a discontinuity in accretion mode [START_REF] Marziani | Searching for the physical drivers of eigenvector 1: influence of black hole mass and Eddington ratio[END_REF][START_REF] Marziani | Lowand high-z highly accreting quasars in the 4D Eigenvector 1 context[END_REF], as the most relevant physical parameter governing the MS and the differences summarized in Table 1 is likely Eddington ratio [START_REF] Marziani | Searching for the Physical Drivers of the Eigenvector 1 Correlation Space[END_REF][START_REF] Shen | The diversity of quasars unified by accretion and orientation[END_REF].

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, [START_REF] Strateva | Double-peaked Low-Ionization Emission Lines in Active Galactic Nuclei[END_REF], redward asymmetric, Fe II opt is weak, the spectral energy distribution (SED) is hard [START_REF] Laor | The Soft X-Ray Properties of a Complete Sample of Optically Selected Quasars[END_REF][START_REF] Shang | Spectral Properties from Lyα to Hα for an Essentially Complete Sample of Quasars[END_REF], [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 [START_REF] Sulentic | C IV λ1549 as an Eigenvector 1 Parameter for Active Galactic Nuclei[END_REF][START_REF] Richards | Unification of Luminous Type 1 Quasars through C IV Emission[END_REF][START_REF] Coatman | C iv emission-line properties and systematic trends in quasar black hole mass estimates[END_REF].

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 [START_REF] Sulentic | The Demise of the Classical Broad-Line Region in the Luminous Quasar PG 1416-129[END_REF] 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).

METHOD

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 [START_REF] Fraix-Burnet | Multivariate approaches to classification in extragalactic astronomy[END_REF] for a review on unsupervised classification techniques in extragalactic astronomy, and to Fraix-Burnet et al. (2006a,b, 2012); [START_REF] Fraix-Burnet | Statistics for Astrophysics: Clustering and Classification[END_REF] 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 [START_REF] Swofford | Paup*: Phylogenetic analysis using parsimony (*and other methods)[END_REF], with a ratchet method to avoid as much as possible local minima [START_REF] Nixon | The parsimony ratchet, a new method for rapid parsimony analysis[END_REF].

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 [START_REF] Boyle | The 2dF QSO Redshift Survey -I. The optical luminosity function of quasi-stellar objects[END_REF] 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., [START_REF] Zamfir | Detailed characterization of Hβ emission line profile in low-z SDSS quasars[END_REF][START_REF] Shen | The diversity of quasars unified by accretion and orientation[END_REF].

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 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, in preparation). The black hole mass (M BH ) has been estimated from Hβ FWHM using the [START_REF] Vestergaard | Determining Central Black Hole Masses in Distant Active Galaxies and Quasars. II. Improved Optical and UV Scaling Relationships[END_REF] 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 [START_REF] Elvis | Atlas of quasar energy distributions[END_REF][START_REF] Richards | Spectral Energy Distributions and Multiwavelength Selection of Type 1 Quasars[END_REF].

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 [START_REF] Fraix-Burnet | The Fundamental Plane of Early-Type Galaxies as a Confounding Correlation[END_REF] 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 M215 sample finally includes seven parameters for the cladistic analysis: radio loudness parameter R K (coded as logRK in the Figures) following [START_REF] Kellermann | VLA observations of objects in the Palomar Bright Quasar Survey[END_REF], R FeII (RFE), FWHM(Hβ) (FWHMHb), the line centroid displacement of Hβ at quarter maximum c( 14 ) (c1o4Hb), W([OIII]λ5007) (WOIII, in Å), the peak shift of [OIII]λ5007 (voIII), the decimal logarithms of bolometric luminosity L bol (logLbol), M BH (logMBH). All radial velocity displacements are referred to the quasar rest frame, in units of km s -1 . Note that voIII has 28 unknown values.

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( 12 ) (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 [START_REF] Hawking | Black hole explosions?[END_REF] 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 and3. 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. 12 ), 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 and5). 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) and CDs (19). Group 18 includes one core-bright LDs and three radio intermediate sources, 0 log R K 1.8. Radio-loud groups are from 39 to 41 in M215 (41 predominantly made of CDs, 40 of LDs). The intermediate groups include low-right corner Pop. A sources but are mainly Pop. B. The cladograms (Figs. 2 and3) trace the relations among the groups. Even if no evolutionary inference should be inferred at this point, it is interesting to note that there is a sequence going from extreme Pop. A and Pop. B RQ, to Pop. B RL. The bottom groups are core-dominated and lobe-dominated RL sources, which are monophyletic groups in M85. In M85 and M215, the different groups appear to follow defined systematic trends, with FWHM(Hβ), R FeII and Γ S decreasing from the first group to the last. Large CIVλ1549 blueshifts seem to be restricted to the first group. The cladistic analysis therefore recovers the trends that define the quasars MS in 4DE1 parameter space involving FWHM(Hβ), R FeII , c( 12 ) CIVλ1549, and Γ S . High values of the L/L Edd are restricted to the first groups with the smallest masses. 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 [START_REF] Urry | Unified Schemes for Radio-Loud Active Galactic Nuclei[END_REF] 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 [START_REF] Wills | Relativistic beaming and quasar emission lines[END_REF][START_REF] Zamfir | New insights on the QSO radio-loud/radio-quiet dichotomy: SDSS spectra in the context of the 4D eigenvector1 parameter space[END_REF][START_REF] Buttiglione | An optical spectroscopic survey of the 3CR sample of radio galaxies with z 0.3 . II. Spectroscopic classes and accretion modes in radio-loud AGN[END_REF]. 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.

DISCUSSION

L/L Edd

Several observational parameters are strongly correlated with L/L Edd [START_REF] Marziani | Searching for the Physical Drivers of the Eigenvector 1 Correlation Space[END_REF][START_REF] Kuraszkiewicz | Emission Line Properties of Active Galactic Nuclei from a Post-COSTAR Hubble Space Telescope Faint Object Spectrograph Spectral Atlas[END_REF][START_REF] Shen | The diversity of quasars unified by accretion and orientation[END_REF]. L/L Edd expresses the ratio between radiation and gravitational forces [START_REF] Netzer | The Effect of Radiation Pressure on Emission-line Profiles and Black Hole Mass Determination in Active Galactic Nuclei[END_REF][START_REF] Marziani | Broadline region physical conditions along the quasar eigenvector 1 sequence[END_REF]. The transition between Pop. A and B occurs at L/L Edd ≈ 0.2 ± 0.1 [START_REF] Marziani | Searching for the physical drivers of eigenvector 1: influence of black hole mass and Eddington ratio[END_REF], 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 [START_REF] Sulentic | C IV λ1549 as an Eigenvector 1 Parameter for Active Galactic Nuclei[END_REF][START_REF] Richards | Unification of Luminous Type 1 Quasars through C IV Emission[END_REF][START_REF] Coatman | C iv emission-line properties and systematic trends in quasar black hole mass estimates[END_REF]. 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).

Ontogenesis and phylogenesis of quasars

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 [START_REF] Mineshige | Slim-Disk Model for Soft X-Ray Excess and Variability of Narrow-Line Seyfert 1 Galaxies[END_REF], 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., [START_REF] Fontanot | The many manifestations of downsizing: hierarchical galaxy formation models confront observations[END_REF][START_REF] Reviglio | Active Galaxies in the Sloan Digital Sky Survey II: galaxy and activity evolution[END_REF][START_REF] Hirschmann | Cosmological simulations of black hole growth: AGN luminosities and downsizing[END_REF]. The very massive black holes that were shining bright mostly belong now to spent systems [START_REF] Lynden-Bell | Galactic Nuclei as Collapsed Old Quasars[END_REF], 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, [START_REF] Walsh | The M87 Black Hole Mass from Gas-dynamical Models of Space Telescope Imaging Spectrograph Observations[END_REF]) 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 [START_REF] Shakura | Black holes in binary systems. Observational appearance[END_REF]. 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 ∆t = c 2 (M BH -M 0 BH ) L η 1 -η ≈ 5.7 • 10 7 M 0 BH,8 (f -1)L -1 46 η 0.1 yr (1)
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 [START_REF] Boyle | The 2dF QSO Redshift Survey -I. The optical luminosity function of quasi-stellar objects[END_REF] luminosity function parameterization:

Φ(M B , z) = Φ(M * B ) 10 0.4[(a+1)(M B -M * B (z))] + 10 0.4[(b+1)(M B -M * B (z))] , (2)
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 [START_REF] Sulentic | VLT/ISAAC spectra of the Hβ region in intermediate-redshift quasars. II. Black hole mass and Eddington ratio[END_REF][START_REF] Natarajan | Is there an upper limit to black hole masses?[END_REF][START_REF] King | How big can a black hole grow?[END_REF], 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:

N lum = z 2 z 1 ñdV = z 2 z 1 M B,max -∞ Φ(M B , z)dM B 4πd 2 C dz E(z) , (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 [START_REF] Zamfir | Detailed characterization of Hβ emission line profile in low-z SDSS quasars[END_REF], 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 [START_REF] Sulentic | GTC spectra of z ≈ 2.3 quasars: comparison with local luminosity analogs[END_REF]. 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 [START_REF] Kelly | Constraints on Black Hole Growth, Quasar Lifetimes, and Eddington Ratio Distributions from the SDSS Broad-line Quasar Black Hole Mass Function[END_REF]. 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). [START_REF] Mathur | Narrow-line Seyfert 1 galaxies and the evolution of galaxies and active galaxies[END_REF] 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 [START_REF] Sulentic | GTC spectra of z ≈ 2.3 quasars: comparison with local luminosity analogs[END_REF]. 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, [START_REF] Wisotzki | The Hamburg/ESO survey for bright QSOs. III. A large flux-limited sample of QSOs[END_REF], objects whose luminosity may also be enhanced by anisotropic emission [START_REF] Urry | The effect of anisotropic emission from thick accretion disks on the luminosity functions of active galactic nuclei[END_REF][START_REF] Dipompeo | The Intrinsic Quasar Luminosity Function: Accounting for Accretion Disk Anisotropy[END_REF]. 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., [START_REF] Wagner | Ultrafast Outflows: Galaxy-scale Active Galactic Nucleus Feedback[END_REF][START_REF] King | Powerful Outflows and Feedback from Active Galactic Nuclei[END_REF], 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 [START_REF] Kaastra | A fast and long-lived outflow from the supermassive black hole in NGC 5548[END_REF], 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., [START_REF] Cavaliere | The Fall of the Quasar Population[END_REF][START_REF] Hopkins | A Unified, Merger-driven Model of the Origin of Starbursts, Quasars, the Cosmic X-Ray Background, Supermassive Black Holes, and Galaxy Spheroids[END_REF]. 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 [START_REF] Zamfir | New insights on the QSO radio-loud/radio-quiet dichotomy: SDSS spectra in the context of the 4D eigenvector1 parameter space[END_REF] suggested a value that exceeds by an order of magnitude the canonical limit introduced by [START_REF] Kellermann | VLA observations of objects in the Palomar Bright Quasar Survey[END_REF]: 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 [START_REF] Bonzini | Star formation properties of sub-mJy radio sources[END_REF][START_REF] Padovani | The faint radio sky: radio astronomy becomes mainstream[END_REF]. 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., [START_REF] Komossa | Radio-loud Narrow-Line Type 1 Quasars[END_REF][START_REF] Berton | Compact steep-spectrum sources as the parent population of flat-spectrum radio-loud narrow-line Seyfert 1 galaxies[END_REF].

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 [START_REF] Marziani | Searching for the physical drivers of eigenvector 1: influence of black hole mass and Eddington ratio[END_REF], 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 [START_REF] Richards | Unification of Luminous Type 1 Quasars through C IV Emission[END_REF][START_REF] Sulentic | C IV λ1549 as an Eigenvector 1 Parameter for Active Galactic Nuclei[END_REF], 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 [START_REF] Sikora | Radio Loudness of Active Galactic Nuclei: Observational Facts and Theoretical Implications[END_REF]. 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 [START_REF] Iwasawa | The variable iron K emission line in MCG-6-30-15[END_REF]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 [START_REF] Blandford | Electromagnetic extraction of energy from Kerr black holes[END_REF], the jet power is related to M BH , the spin angular momentum J, and the magnetic field by:

P ν = 5 • 10 10 J J max 2 M BH 1M 2 B 1G 2 erg s -1 Hz -1 , (4)
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 [START_REF] Punsly | Evidence of the Dynamics of Relativistic Jet Launching in Quasars[END_REF] 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 [START_REF] Narayan | Magnetically Arrested Disk: an Energetically Efficient Accretion Flow[END_REF]. 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., [START_REF] Blandford | Hydromagnetic flows from accretion discs and the production of radio jets[END_REF][START_REF] Riols | Magnetorotationally driven wind cycles in local disc models[END_REF].

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 [START_REF] Wu | The black hole mass, Eddington ratio and M BH -σ [OIII] relation in young radio galaxies[END_REF] 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 , Sulentic et al. 2015), is very powerful at present and also shows evidence of relic radio emission, probably due to a past activity cycle [START_REF] Punsly | The extreme ultraviolet spectra of low-redshift radio-loud quasars[END_REF]. 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 [START_REF] Boroson | The emission-line properties of low-redshift quasi-stellar objects[END_REF] 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. 1: [START_REF] Boroson | The emission-line properties of low-redshift quasi-stellar objects[END_REF]2: Sulentic et al. 2000a;3: Collin et al. 2006;4: Shen and Ho 2014;5: Sulentic et al. 2007;6: Baskin and Laor 2005;7: Richards et al. 2011;8: Sulentic et al. 2016;9: Wang et al. 1996: 10: Bensch et al. 2015;11: Véron-Cetty et al. 2001;12: Sulentic et al. 2002;13: Marziani et al. 2003b;14: Marziani et al. 2001;15: Wills et al. 1999;16: Bachev et al. 2004;17: Coatman et al. 2016;18:Zhang et al. 2011;19: Marziani et al. 2016a;: Zamanov et al. 2002;21: Wang et al. 2006;22: Turner et al. 1999;23: Grupe et al. 2001;24: Giveon et al. 1999;25: Zamfir et al. 2008;26: Reichard et al. 2003;27: Sulentic et al. 2006;28: Negrete et al. 2012;29: Boroson 2002;30: Peterson et al. 2004;31: Kuraszkiewicz et al. 2000.

  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. (

Figure 1 .

 1 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).

Figure 2 .

 2 Figure 2. Cladograms for the M85 sample rooted with the lowest-M BH group, showing barplots for the Eddington ratio L/L Edd in addition to the eleven parameters used to find the tree. Colors correspond to the groups defined from the tree on the left.

Figure 3 .

 3 Figure 3. Cladogram for the M215 sample rooted with the lowest-M BH group, showing barplots for L/L Edd , Γ S , CIVλ1549, c(12 ), 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

Figure 4 .

 4 Figure 4. Box plots with identified groups for the M85 sample. The ordinate axis shows the parameter coded on top of the plot, with code meaning explained in §2.3. Medians are thick horizontal lines while boxes define limits from first to third interquartile range. Outliers are shown as open circles. Colors correspond to those in Fig. 2.

Figure 5 .

 5 Figure 5. Box plots with identified groups for the M215 sample. Colors correspond to those in Fig. 3, and the ordinate axis meaning is the same of Fig. 4.

Figure 6 .

 6 Figure 6. Unrooted cladistic trees for the M85 (left) and M215 (right) samples. The main groups are identified by shaded areas in the M85 diagram. Colors corresponds to those in Fig. 2 (M85) and Fig. 3 (M215).

Table 1 .

 1 MAIN TRENDS ALONG THE 4DE1 SEQUENCE

	Parameter	Population A	Population B	References
	FWHM(Hβ BC )	800 -4000 km s -1	4000 -10000 km s -1	1, 2, 3, 4
	R Fe	0.7	0.3	1, 2
	c( 1 2 ) CIVλ1549 BC	-800 km s -1	-250/+70 (RQ/RL)	5, 6, 7, 8
	Γ S	often large (> 2)	rarely large (≈ 2)	2, 9, 4, 10
	W(Hβ BC )	∼ 80 Å	∼ 100 Å	2
	Hβ BC profile shape	Lorentzian	double Gaussian	11, 12, 13
	c( 1 2 ) Hβ BC	∼ zero	+500 km s -1	13
	SiIII / CIII]	0.4	0.2	14, 15,16
	FWHMCIVλ1549 BC	(2-6) •10 3 km s -1	(2-10) •10 3 km s -1	5, 17
	W(CIVλ1549 BC )	few Å -≈ 60 Å	∼ 100 Å	4, 6, 7
	AI(CIVλ1549 BC )	-0.1	0.05	5
	W([OIIIλ5007)	1 -20	20 -80	1, 18, 19
	vr([OIIIλ5007)	negative / 0	∼ 0 km s -1	4, 18, 19, 20
	FIR color α(60, 25)	0 --1	-1 --2	21
	X-ray variability	extreme/rapid common	less common	22, 23
	optical variability	possible	more frequent/higher amplitude 24
	probability radio loud	≈ 3-4%	≈ 0.25 %	4, 25
	BALs	extreme BALs	less extreme BALs	36,37
	log density 1	11	9.5	14, 28
	log U 1	-2.0/-1.5	-1.0/-0.5	14, 28
	log M BH [M ]	6.5 -8.5	8.0 -9,5	7, 8, 29
	L/L Edd	≈ 0.2 -1.0	∼ 0.01 -≈ 0.2	1, 4, 7, 29, 30, 31

This is a provisional file, not the final typeset article

Frontiers

http://cdsarc.u-strasbg.fr/viz-bin/Cat?J/ApJS/145/199

http://cdsarc.u-strasbg.fr/viz-bin/Cat?J/ApJ/666/757

The group identification is reported for each quasar in the file quasargroups.csv appended to the paper, Table2.Frontiers

Refers to a dense, low ionization emitting region that is apparently varying along the E1 sequence, and is associated with the reverberation response of Hβ(Negrete et al.

2013).

ACKNOWLEDGMENTS

FUNDING

DD acknowledges support from grants PAPIIT108716, UNAM, and CONACyT221398.

CONFLICT OF INTEREST STATEMENT

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Laurent U Perrinet
email: laurent.perrinet@univ-amu.fr

Biologically-inspired characterization of sparseness in natural images

Keywords: Machine vision, Image texture, Computer vision, Neuroscience

Natural images follow statistics inherited by the structure of our physical (visual) environment. In particular, a prominent facet of this structure is that images can be described by a relatively sparse number of features. We designed a sparse coding algorithm biologically-inspired by the architecture of the primary visual cortex. We show here that coefficients of this representation exhibit a heavy-tailed distribution. For each image, the parameters of this distribution characterize sparseness and vary from image to image. To investigate the role of this sparseness, we designed a new class of random textured stimuli with a controlled sparseness value inspired by our measurements on natural images. Then, we provide with a method to synthesize random textures images with a given statistics for sparseness that matches that of some given class of natural images and provide perspectives for their use in neurophysiology.

author = {Perrinet, Laurent U.}, booktitle = {2016 6th European Workshop on Visual Information Processing (EUVIP)}, doi = {10.1109/EUVIP.2016.7764592}, isbn = {978-1-5090-2781-1}, month = {oct}, pages = {1--6}, publisher = {IEEE}, title = {Biologically-inspired characterization of sparseness in natural images}, url = {http://ieeexplore.ieee.org/document/7764592/}, year = {2016} }

Motivation

Figure 1: Sparse coding of natural images using log-Gabor edges and applied to "Lena": An instance of the reconstruction of a natural image as a sum of elementary edges using sparse coding. Parameters for each edge are its position, orientation, scale and scalar amplitude of the coefficient. Edges are overlaid on this reconstructed image as colored segments, length denoting scale and hue the orientation -see SI Section 5 for a full description of the algorithm. Edges outside the dashed circle are discarded to avoid artifacts. Few edges (here 1024) are necessary to efficiently represent the image: the distribution of features is sparse. The position and orientation of the set of edges gives information about the shape of visual objects, while the distribution of coefficients characterizes the texture content of the image.

Natural images, that is, visual scenes that are relevant for an animal, most generally consist of the composition of visual objects, from vast textured backgrounds to single isolated items. By the nature of the structure of the physical (visual) world, these visual items are often sparsely distributed, such that these elements are clustered and a large portion of the space is void: A typical image thus contains vast areas which are containing little information while most information is concentrated in a small portion of the image (see Figure 1 and Supplementary Figure 5 ). A first formalism to quantify this sparseness is to consider the distribution of coefficients obtained by coding such images using a bank of filters resembling the edge-like profiles observed in the primary visual cortex (V1) of primates (See Supplementary Figure 4). From this linear representation, it is possible to determine a near-to-optimal coding formalism which involves the inversion of a model of the (forward) neural transformation [START_REF] Simoncelli | Characterization of neural responses with stochastic stimuli[END_REF]. As a matter of fact, a normative explanation for the coding of images into activity in the neural tissue is the optimal representation of visual information [START_REF] Atick | Could information theory provide an ecological theory of sensory processing?[END_REF]. Inspired by the representation in V1, a popular method requires the efficient modeling of images as the sparsest combination of edges and such a coding strategy indeed resembles the neural activity (that is, the spike patterns) observed in V1 [START_REF] Perrinet | Sparse models for computer vision[END_REF]. Using such a method, it was previously shown that any image from a database of static, grayscale natural images may be coded using solely a relatively few number of coefficients (See Figure 2 and SI Section 5 for a full description of the algorithm). Crucially, applying this sparse coding method on this database, we observed that some images were sparser than others but that globally they all fitted well a similar prototypical probability distribution function (pdf). The parameters of this pdf provide with a quantitative descriptor for sparseness and is characteristic for this property. However, it is largely unclear how this property inherent to the structure of natural images may be used in the early visual system or in computer vision.

The different levels of sparseness observed in images are qualitatively linked to their textural content. Indeed, different levels of roughness in any image correspond to different distributions of the coefficients corresponding to its representation. Such an analysis was previously performed for instance by characterizing the Hurst parameter for two-dimensional fractional Brownian motion [START_REF] Lance M Kaplan | Extended fractal analysis for texture classification and segmentation[END_REF]. However, this analysis was done on the linear coefficients of Gabor features, while we rather wish here to characterize sparseness in the above-mentioned biologically inspired generative model for natural images. Compared to a linear representation, this coding is non-linear and the distribution of edges' parameters are a priori different from that described in [START_REF] Lance M Kaplan | Extended fractal analysis for texture classification and segmentation[END_REF]. In particular, while the distribution of position, scale and orientation are well characterized, an unified theory for the distribution of sparse coefficients is still lacking. As such, a major challenge is to refine the characterization of sparseness in natural images such that it conforms to the widest variety of visual scenes. Here, by using a bio-inspired representation based on the architecture of the primary visual cortex, we will first characterize the pdf of sparse coefficients on an image-by-image basis and then propose a simple characterization of sparseness based on the number of non-zero coefficients that we apply to texture synthesis.

Our first contribution is to show that while the distribution of coefficients in natural images is stereotyped, the underlying parameters display some variability (see Figure 2). Our second contribution is an efficient extension of a texture synthesis model which allows to parameterize the sparseness in the texture. This formulation is important to challenge the effects of the model's parameters in future neurophysiological experiments. Indeed, in neurophysiology, a popular technique to challenge the coding hypothesis is to test the system (V1) using complex stimuli in the form of stochastic stimuli [START_REF] Touryan | Analysis of sensory coding with complex stimuli[END_REF]. These stimuli correspond to instances of the generative model for natural images which most likely describe a wide range of image instances. From that perspective, we motivate the generation of an optimal stimulation within a stationary Gaussian dynamic texture model (see Figure 3). We base our model on a previously defined heuristic coined "Motion Clouds" [START_REF] Sanz-Leon | Motion clouds: model-based stimulus synthesis of natural-like random textures for the study of motion perception[END_REF]. These contributions show that overall, sparseness in natural images is an important parameter to understand natural scenes. Beyond the contribution to the understanding of the models underlying image representation which may have applications to tune computer vision algorithms, this may have applications outside computer vision, for instance for the optimal stimulation of neurophysiological settings as for instance in a retinal implant. 1.5 2.

# A B

Figure 2: Coefficients from the Sparse coding of images in the primary visual cortex (V1) follows regular statistics: (A-inset) Coefficients decrease as a function of their rank (average across images ± standard deviation). This decrease is faster when using a sparse coding mechanism (Matching Pursuit, 'MP') than when just ordering the linear coefficients ('lin'). (A) For the sparse coding, we controlled the quality of the reconstruction from the edge information such that the residual energy is less than 5% over the whole set of images, a criterion met on average when identifying 2048 edges per image for images of size 256 × 256 (that is, a relative sparseness of ≈ 0.02% of activated coefficients). (B) When observing on an image-by-image basis the histogram of the coefficients' amplitude, each follows a generic log-Normal probability density function, only controlled by two parameters (the mode and bandwidth B) for each of natural image. While the mode parameter did not vary significantly, higher values of B qualitatively correspond to sparser structures.(B-inset) The distribution of these sparseness parameters shows a variety of values from the most sparsely distributed (right, B ≈ 1.) to those containing mostly dense textures (left, B ≈ 3.).

Biologically-inspired sparse coding

Let us first describe the method used for efficiently extracting edges in natural images before measuring the statistics of the edges' parameters. The first step of our method involves defining the dictionary of templates (or filters) for detecting edges. In particular, we use a synthesis model for the edge representation, so that the edges we detect are guaranteed to be sufficient to regenerate the image with a low error. For that, we use a log-Gabor representation, which is well suited to represent a wide range of natural images [START_REF] Fischer | Sparse approximation of images inspired from the functional architecture of the primary visual areas[END_REF]. This representation gives a generic model of edges parameterized by their shape, orientation, and scale. We set the range of these parameters to match what has been reported for the responses in primates' V1. In particular, we set the bandwidth of the Fourier representation of the filters to 1 and π/8 respectively in log-frequency and polar coordinates to get a family of elongated and thus orientation-selective filters (see [START_REF] Fischer | Self-invertible 2D log-Gabor wavelets[END_REF] , see Supplementary Figure 4). Prior to the analysis of each image, we used the spectral whitening filter [START_REF] Olshausen | Sparse coding with an overcomplete basis set: A strategy employed by V1?[END_REF] to provide a good balance of the energy of output coefficients [START_REF] Fischer | Sparse approximation of images inspired from the functional architecture of the primary visual areas[END_REF][START_REF] Perrinet | Coding static natural images using spiking event times: Do neurons cooperate[END_REF]. Although orthogonal representations are popular for computer vision due to their computational tractability, it is desirable in our context that we have a high over-completeness in the representation. Indeed, the representation of edges, similar to the receptive fields in the primary visual cortex, should provide with a relative invariance to common geometrical transformations (translation, rotation, and scale). First, a linear convolution model automatically provides a translation-invariant representation. Such invariance can be extended to rotations and scalings by choosing to multiplex these sets of filters at different orientations and spatial scales. Ideally, the parameters of edges would vary in a continuous fashion, to provide relative invariance. For the set of 256 × 256 images, we chose to have 8 dyadic levels (that is, doubling the scale at each level) with 24 different orientations. Orientations are measured as an undirected angle in radians, in the range from 0 to π (but not including π). Tests with a range of different numbers of orientations and scales yielded similar results [START_REF] Perrinet | Sparse models for computer vision[END_REF]. Finally, each image is transformed into a pyramid of coefficients. This pyramid consists of approximately 4/3 × 2562 ≈ 8.7 × 10 4 pixels multiplexed on 8 scales and 24 orientations, that is, approximately 16.7 × 10 6 coefficients, an over-completeness factor of about 256.

This transform is linear and can be performed by a simple convolution repeated for every edge type. Following [START_REF] Fischer | Self-invertible 2D log-Gabor wavelets[END_REF], convolutions were performed in the Fourier (frequency) domain for computational efficiency. The Fourier transform allows for a convenient definition of the edge filter characteristics, and convolution in the spatial domain is equivalent to a simple multiplication in the frequency domain. By multiplying the envelope of the filter and the Fourier transform of the image, one may obtain a filtered spectral image that may be converted to a filtered spatial image using the inverse Fourier transform. We exploited the fact that by omitting the symmetrical lobe of the envelope of the filter in the frequency domain, the output of this procedure gives a complex number whose real part corresponds to the response to the symmetrical part of the edge, while the imaginary part corresponds to the asymmetrical part of the edge (see [START_REF] Fischer | Self-invertible 2D log-Gabor wavelets[END_REF] for more details). More generally, the modulus of this complex number gives the energy response to the edge (comparable to the response of complex cells in area V1), while its argument gives the exact phase. This property further expands the richness of the representation. Overall, such a representation is implemented in the LogGabor package of Python scripts 1 .

Because this dictionary of edge filters is over-complete, the linear representation would give an inefficient representation of the distribution of edges (and thus of edge co-occurrences) due to a priori correlations between coefficients. Therefore, starting from this linear representation, we searched for the most sparse representation. Minimizing the 0 pseudo-norm (the number of non-zero coefficients) leads to an expensive combinatorial search with regard to the dimension of the dictionary (it is NP-hard). As proposed first by [START_REF] Perrinet | Sparse spike coding in an asynchronous feed-forward multi-layer neural network using matching pursuit[END_REF], we may approximate a solution to this problem using a greedy approach, namely the Matching Pursuit algorithm. This class of algorithms gives a generic and efficient representation of edges, as illustrated by the example in Figure 1. For the simulations presented here, sparse coding of images is implemented in the SparseEdges package of Python scripts 2 which depend on packages NumPy (version 1.11.1) and SciPy (version 0.18.0) [START_REF] Oliphant | Python for scientific computing[END_REF] on a cluster of Linux computing nodes. Visualization was performed using Matplotlib (version 1.5.1) [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF]. Our goal is to study how the statistics of edges' occurrences vary across a set of natural images. In order to perform that, wes first defined a database of 600 images 3 . The performance of the algorithm can be measured quantitatively by reconstructing the image from the list of extracted edges. Measuring the ratio of extracted energy in the images, N = 2048 edges were enough to extract on average more than 95% of the energy of 256 × 256 images on all sets of images (see Figure 2-A). Then, we tested different hypothesis for the distribution of coefficients and tested the synthesis of textures using this generative model well tuned to these natural images. All these different steps are reproducible as a set of notebooks.

3 Distribution of sparse coefficients and the design of the "DropLets" stimuli Empirically, we observed that the distribution of sparse coefficients for any natural image fits with a log-Normal probability distribution function (see Figure 2-B). In particular, the histogram of occurrences of the absolute sparse coefficients for any given image follows a smooth curve in a log-log plot, suggesting a fit with a log-Nornal distribution as parameterized by the mode and bandwidth B. We fitted each histogram using the Levenberg-Marquardt method from SciPy [START_REF] Oliphant | Python for scientific computing[END_REF].

In particular, this was validated by checking the goodness of fit for each image. While the mode did not vary significantly across the database of natural images, we observed that the bandwidth B of the distribution was significantly different across the database. This is shown as the histogram of the B values obtained for each image (see Figure 2-B, inset). This demonstrates that natural images show a wide range of sparseness in their coefficients from very sparse (B ≈ 3.) to more dense (B ≈ 1.). Qualitatively, this corresponds to images which are containing respectively less to more dense textures. An interpretation for that class of distributions comes from the fact that such pdf emerges from the mixing of different sources with narrow bandwidths: a denser mixing produces a broader bandwidth. As a consequence, sparseness could be simply parameterized as the intensity of the mixing of such events, or put even more simply, to the number of non-zero coefficients. Such knowledge can now be injected in a model for texture synthesis. From these observations, we propose a mathematically-sound derivation of a general parametric model of synthetic textures. This model is defined by aggregation, through summation, of a basic spatial "texton" template φ(x, y). The summation reflects a transparency hypothesis, which has been adopted for instance in [START_REF] Galerne | Micro-Texture synthesis by phase randomization[END_REF]. In particular, this simple generative model considers independent, transparent elementary features. While one could argue that this hypothesis is overly simplistic and does not model occlusions, it leads to a tractable framework of stationary Gaussian textures, which has proved useful to model static micro-textures [START_REF] Galerne | Micro-Texture synthesis by phase randomization[END_REF] and dynamic natural phenomena [START_REF] Xia | Synthesizing and mixing stationary gaussian texture models[END_REF]. In particular, the simplicity of this framework allows for a fine tuning of frequency-based (Fourier) parameterization, which is desirable for the interpretation of neuro-physiological experiments. In summary, it states that luminance I(x) for x ∈ R 2 is defined = 5e 06 = 1e 04 = 2e 03 = 5e 02 = 1e + 00 where, similarly to the variety observed in a set of natural images, the distribution of coefficients is parameterized by different levels of sparseness. From the observations in a range of natural images (see Figure 2-B, inset), we define the sparseness as the relative number of non-zero coefficients, from very sparse ( = 5e -06) to dense ( = 1), and with values chosen on a geometrical scale between those extremes. Note that, while for = 1 these textures are fully dense and correspond to a linearly filtered noise, other values correspond to progressively sparser textures as tends to lower values. Moreover, while the textons used to synthesize these textures are arbitrary in our formalism, targetting a neurophysiological experiment on the retina, we have chosen here for the sake of simplicity to use symmetrical Mexican-hat profiles ("drops" , see Supplementary Figure 4) at a single scale which is optimal for retinal ganglion cells.

as a random field [START_REF] Vacher | Biologically inspired dynamic textures for probing motion perception[END_REF]:

I(x) = i∈N a i • φ(x -x i ) (1)
where the respective random variables parameterize respectively each texton's scalar value a i and position x i ∈ R 2 . As was previously mentioned [START_REF] Sanz-Leon | Motion clouds: model-based stimulus synthesis of natural-like random textures for the study of motion perception[END_REF][START_REF] Vacher | Biologically inspired dynamic textures for probing motion perception[END_REF], this defines the set of "Motion Clouds" stimuli as that defined by:

I(x) = [ i∈N a i • δ(x -x i )] * φ(x) (2)
where * denotes the convolution operator over variable x. Noting the Fourier transform as F and f ∈ R 2 the frequency variables, the image I is thus a stationary Gaussian random field of covariance having the power-spectrum E = F( φ) (for a proof, see [START_REF] Vacher | Biologically inspired dynamic textures for probing motion perception[END_REF]). It comes

F(I)(f ) = A • e i•Φ • E(f ) where A • e i•Φ = i∈N a i • e -2iπ f, xi
corresponds to an iid random phase field scaled by A ∈ R + . Such a random texture is parameterized by the (positive-, real-valued) envelope E. To match the statistics of natural scenes or some category of textures, the envelope E is usually defined as some average spatio-temporal coupling over a set of natural images. Note that it is possible to consider any arbitrary texton F -1 (E), which would give rise to more complicated parameterizations for the power spectrum.

In particular, we may use the same set of edge-like filters used in the sparse coding algorithm above.

The originality provided here is to introduce an explicit sparsity in equation 2 by parameterizing a sparse distribution for the coefficients a i . Mathematically, the a i are drawn from a heavy-tailed probability distribution function p a such as that found in natural images (See Figure 2-B). Such assumptions were previously used for generating procedural noise in computer vision [START_REF] Ares Lagae | Procedural noise using sparse gabor convolution[END_REF], but we focus here on a definition based on a model of sparseness in images. A similar endeavor was initiated by [START_REF] Sallee | Learning Sparse Multiscale Image Representations[END_REF] by defining a mixture of a Dirac and a Gaussian distributions, but here, we derive it from the sparse coding observed in a set of natural images. First, let's define E = i a i • δ(x -x i ) the (sparse) matrix corresponding to the textons' coefficients at different positions. We also note that we can still compute the image in Fourier space, but with a random phase field which by linearity of F is simply given ∀f by :

F(I)(f ) = F(E)(f ) • E(f ) (3)
Compared to the pure summation, a first advantage of this procedure is that the computation time does not depend on the number of events. A further advantage of equation 3 in generating these stimuli is that for any given instance of the noise, we know the position (x, y) of each event. When using rotationally symmetric, Mexican-hat profiles, the resulting texture resembles the shot noise image of drops of rain on water, such that we coined this set of stimulus as "DropLets" (See Figure 3). In all generality, for any given texton (for instance a choice of a given orientation θ and scale σ), the implementation of this texture generation is of the same order of complexity as that of Motion Clouds and allows for the easy control of the sparseness in the stimulus. These different (uncoherent) "layers" may then be summed to form a texture of sparsely distributed edge-like features.

Perspectives in neurophysiology

In this paper, we have characterized the different levels of sparseness which are present in natural images and synthesized random textures with parameterized levels of sparseness. First, we have replicated the observation that natural images follow a prototypical structure for the probability density function of the coefficients that characterize them. Importantly, we have shown that on an image-by-image basis, this structure is qualitatively well captured by using the number of non-zero coefficients as a measure of the sparseness of a given image, from a dense texture to a highly sparse configuration. Based on these results, we designed random textures which replicate this parameterization of their sparseness. In a recent experimental project using such a set of textures, it has been shown that the neural code as recorded in the retina responds differentially to these different levels of sparseness [START_REF] Ravello | Spatiotemporal tuning of retinal ganglion cells dependent on the context of signal presentation[END_REF]. In particular, we used the set of stimuli (See Figure 3) to evaluate the response on ganglion cells on an ex vivo preparation and analyzed the response of the population of neurons as a function of the sparseness parameter. Such analysis proved that the retinal neural code showed differential responses with respect to the sparseness of images and preliminary measures on efficiency suggest that it is tuned for a given level of sparseness: In that experimental results, we have shown that there is a limit in sparseness for which the retina can respond optimally; beyond that point, the response is more or less the same, meaning that the retina is still coding the features present in the sparser sequences but is not responding to the additional features of the more dense sequences. We can relate these results to the hypothesis of optimal coding of natural scenes [START_REF] Wilson | Visual Perception and the Statistical Properties of Natural Scenes[END_REF], in which the visual system has a limited capacity to transmit information, it has adapted through evolution to a small subset of all the possible images and optimally code them, discarding the irrelevant information. Thus, the efficiency of the retina to code the stimuli reaches its peak at a sparseness level that would be closer to what the system has evolved to code, and when presented with stimuli containing more signals, it discards the additional information. Much research has been performed to investigate the relationship between natural images and optimal coding, although the focus has been mainly on the spatiotemporal correlations [START_REF] Pitkow | Decorrelation and efficient coding by retinal ganglion cells[END_REF][START_REF] Sur | Spatial Correlations in Natural Scenes Modulate Response Reliability in Mouse Visual Cortex[END_REF]. We have observed that when keeping the spatiotemporal components constant, the modulation of the sparseness of the stimuli has an effect on the retinal response, and more importantly, it allows us to see the level of sparseness beyond which the efficiency does not increase.

It is important to note that these results are the outcome of an interdisciplinary convergence between image processing (to characterize sparseness in natural images), mathematical modeling (for the synthesis of textures) and neurophysiology (for the recordings and their analysis). In particular, in this original framework, neural recordings are not analyzed post-hoc, but are instead tuned by the design of parameterized stimuli. In particular, these stimuli are defined from the analysis of natural images. One limit of the neurophysiological study is that we limited ourselves to a simplistic class of textons (Mexican-hat shaped profiles , compare to Supplementary Figure 4), both for the analysis and synthesis but that ganglion cells in the retina are known to be selective to a wider class of stimulations. However, we believe that this class of stimuli is general enough to characterize a wide range of different cell types. Indeed, by looking at local combinations of events (such as doublets), one could characterize different sub-types of ganglion cells, both static, oriented or even moving, such as to characterize for instance directionally selective ganglion cells. In particular, by manipulating the statistics in the event's matrix E, one could target more specifically each of these sub-types. Moreover, while it was found that such textures would be useful in characterizing the response of retinal neurons, this should prove also very useful to characterize the response of V1 neurons. In particular, it would be essential to tune the relation of "neighboring edges" as characterized by how close they are in visual space space but also in orientation and scale, as measured by the so-called "association field" [START_REF] Perrinet | Edge co-occurrences can account for rapid categorization of natural versus animal images[END_REF].

Supplementary material: Sparse coding algorithm

In general, a greedy approach is applied when finding the best combination is difficult to solve globally, but can be solved progressively, one element at a time. Applied to our problem, the greedy approach corresponds to first choosing the single filter Φ i that best fits the image along with a suitable coefficient a i , such that the single source a i Φ i is a good match to the image. Examining every filter Φ j , we find the filter Φ i with the maximal correlation coefficient, where:

i = argmax j I I , Φ j Φ j , (4)
•, • represents the inner product, and • represents the 2 (Euclidean) norm. Since filters at a given scale and orientation are generated by a translation, this operation can be efficiently computed using a convolution, but we keep this notation for its generality. The associated coefficient is the scalar projection:

a i = I, Φ i Φ i 2 (5)
Second, knowing this choice, the image can be decomposed as

I = a i Φ i + R (6)
where R is the residual image. We then repeat this 2-step process on the residual (that is, with I ← R) until some stopping criterion is met. Note also that the norm of the filters has no influence in this algorithm on the choice function or on the reconstruction error. For simplicity and without loss of generality, we will thereafter set the norm of the filters to 1: ∀j, Φ j = 1. Globally, this procedure gives us a sequential algorithm for reconstructing the signal using the list of sources (filters with coefficients), which greedily optimizes the 0 pseudo-norm (i.e., achieves a relatively sparse representation given the stopping criterion). The procedure is known as the Matching Pursuit (MP) algorithm [START_REF] Mallat | Matching Pursuit with time-frequency dictionaries[END_REF], which has been shown to generate good approximations for natural images [START_REF] Perrinet | Role of homeostasis in learning sparse representations[END_REF].

For this work we made two minor improvements to this method: First, we took advantage of the response of the filters as complex numbers. As stated above, the modulus gives a response independent of the phase of the filter, and this value was used to estimate the best match of the residual image with the possible dictionary of filters (Matching step). Then, the phase was extracted as the argument of the corresponding coefficient and used to feed back onto the image in the Pursuit step. This modification allows for a phase-independent detection of edges, and therefore for a richer set of configurations, while preserving the precision of the representation.

Second, we used a "smooth" Pursuit step. In the original form of the Matching Pursuit algorithm, the projection of the Matching coefficient is fully removed from the image, which allows for the optimal decrease of the energy of the residual and allows for the quickest convergence of the algorithm with respect to the 0. This is often a problem when edges are aligned (e.g. on a smooth contour), as the different parts will be removed independently, potentially leading to a residual with gaps in the line. Our goal here is not to get the fastest decrease of energy, but rather to provide a good representation of edges along contours. We therefore used a more conservative approach, removing only a fraction (denoted by α) of the energy at each pursuit step (for MP, α = 1). We found that α = 0.5 was a good compromise between rapidity and smoothness. One consequence of using α < 1 is that, when removing energy along contours, edges can overlap; even so, the correlation is invariably reduced. Higher and smaller values of α were also tested, and gave classification results similar to those presented here.

In summary, the whole learning algorithm is given by the following nested loops in pseudo-code:

1. draw a signal I from the database; its energy is E = I 2 , 2. initialize sparse vector s to zero and linear coefficients ∀j, a j =< I, Φ j >,

Figure 3 :

 3 Figure 3: DropLets stimuli: We design DropLets as random phase textureswhere, similarly to the variety observed in a set of natural images, the distribution of coefficients is parameterized by different levels of sparseness. From the observations in a range of natural images (see Figure2-B, inset), we define the sparseness as the relative number of non-zero coefficients, from very sparse ( = 5e -06) to dense ( = 1), and with values chosen on a geometrical scale between those extremes. Note that, while for = 1 these textures are fully dense and correspond to a linearly filtered noise, other values correspond to progressively sparser textures as tends to lower values. Moreover, while the textons used to synthesize these textures are arbitrary in our formalism, targetting a neurophysiological experiment on the retina, we have chosen here for the sake of simplicity to use symmetrical Mexican-hat profiles ("drops" , see Supplementary Figure4) at a single scale which is optimal for retinal ganglion cells.

Figure 4 :

 4 Figure4: Filters used in this paper. To mimic the average profile of receptive fields in the primary visual cortex, we use localized, orientation-selective filters. These are characterized by a cone-shaped spectrum (Left) and provide with the log-Gabor wavelets (Right)

3 .

 3 while the residual energy E = I 2 is above a given threshold do: (a) select the best match: i = ArgMax j |a j |, where | • | denotes the modulus, (b) increment the sparse coefficient:s i = s i + α • a i , (c) update residual image: I ← I -α • a i • Φ i , (d) update residual coefficients: ∀j, a j ← a j -α • a i < Φ i , Φ j >,4. the final non-zero values of the sparse representation vector s, give the list of edges representing the image as the list of couples (i, s i ), where i represents an edge occurrence as represented by its position, orientation and scale.

Figure 5 :

 5 Figure 5: Progressive reconstruction as a function of sparseness.

Available at https://github.com/bicv/LogGabor and documented at https://pythonhosted.org/LogGabor.

These scripts are available at https://github.com/bicv/SparseEdges and documented at https://pythonhosted.org/SparseEdges.

Publicly available at http://cbcl.mit.edu/software-datasets/serre/SerreOlivaPoggioPNAS07.

pseudo-norm (i.e., it rapidly achieves a sparse reconstruction with low error). However, this efficiency comes at a cost, because the algorithm may result in non-optimal representations due to choosing edges sequentially and not globally.