The situation is more complicated for galaxy models. The monolithic hydrodynamic models by (Chiosi and Carraro, 2002), shortly indicated CC-A and CC-B and the early-hierarchical models by (Merlin et al., 2012), shortly indicated M-M] provide the following MRRs: log R e = 0.331 ⁡ log M s − 3.644  CC - A (5) log R e = 0.273 ⁡ log M s − 1.994  CC - B (6) log R e = 0.241 ⁡ log M s − 1.750  M - M  (7)

We recall that the three groups of models (identical in the input physics) are calculated with different formation redshift z _f (hence initial density): CC-A have z _f ≃ 5, CC-B z _f ≃ 1, and M-M z _f ≃ 1 − 2. In the MR plane they lay on lines with similar slope but different zero points. This suggests that the slope is linked to the physical structure of the models while the zero-point is reminiscent of the initial density. Surprisingly, the slopes of the above relations are not identical to that of ETGs (Eq. 4), but close to that of DGs. Furthermore, along the sequence of each group, the duration of the star formation activity is long and in a burst-like mode of low intensity in low mass galaxies and short and intense (often a single burst of activity) in the high mass ones. Remarkably, only the most massive galaxies formed in redshift interval 5 ≥ z _f ≥ 2, in which star formation has ceased long ago, may fall into the region of ETGs.

The Illustris hierarchical models provide similar relationships, once they are split into two groups: log R e = 0.297 ⁡ log M s − 2.513 f o r log M s ≤ 10.5 , (8) log R e = 0.519 ⁡ log M s − 4.492 f o r log M s ≥ 10.5 . (9)

The first relation holds for the vast majority of models and reminds that of normal DGs, while the second one holds for a small group of objects and is close to the case of ETGs. In the hierarchical scheme the models of the first group (in Eq. 8) are the seeds of those in the second group located along the MRR of Eq. 9.

Finally, there is the MRR proposed by Fan et al. (2010). This is derived in the following way. Independently of the monolithic or hierarchical scheme, the seeds of galaxies are perturbations of matter made of DM and BM. These collapse when the density contrast with respect to the surrounding medium reaches a suitable value. Assuming spherical symmetry and indicating with M _T and R _T the total mass and associated radius, and making the approximation M _T = M _D + M _B ≃ M _D and R _T ≃ R _D , the mass-radius relation for each individual galaxy is R D 3 = 3 4 π M D λ ρ u ( z ) → R D ∝ M D 1 / 3 1 + z f (10) where ρ u ( z f ) ∝ ( 1 + z f ) 3 is the density of the Universe at the collapse redshift z _f , and λ the density contrast of the DM halo. This expression has a general validity, whereas λ depends on the cosmological model of the Universe, including the ΛCDM case. All details and demonstration of it can be found in Bryan and Norman [Bryan and Norman (1998), their Eq. 6]. The collapse increases the mean density of DM and BM so that, when a critical value of the BM density is reached, stars can form at the center of the system under suitable star formation rates. In the context of the ΛCDM cosmology, Fan et al. (2010) have adapted the general relation (10) to provide an equation connecting the halo mass M _D , the stellar mass M _s , the half light (mass) radius R _e , the shape of the BM S _S (n _S ) related to the Sérsic profile index n _S , the velocity dispersion f _σ of the BM component with respect to that of DM, and finally the ratio m = M _D /M _s . The expression is R e = 0.9 S S ( n S ) 0.34 25 m 1.5 f σ 2 M D 1 0 12 M ⊚ 1 / 3 4 ( 1 + z f ) . (11) where f _σ yields the three dimensional stellar velocity dispersion as a function of the DM velocity dispersion σ _s = f _σ σ _D (here we adopt f _σ = 1). The typical value for S _S (n _S ) is 0.34. For more details see Fan et al. (2010) and references therein.

The most important parameter of Eq. 11 is the ratio m = M _D /M _s . Using the Illustris data Chiosi et al. (2019) investigated how this ratio varies in the mass interval 8.5 < log  M _D < 13.5 (masses are in M _⊙) and from z = 0 to z = 4 (see Section 11). They find that the following relation is good for all practical purposes log ⁡ m = log M D M s = 0.062 log M D + 0.429 . (12)

The slope of the Fan et al. (Fan et al., 2010) relation, which visualizes the position on the MR plane of systems born at the same redshift once their stars are formed, is 0.333. This is very similar to that of theoretical models, i.e. Eqs. 5–8.

The most intriguing question to answer is “Why is the observational MRR for ETGs so different from the theoretical one?”

The observed distribution of astrophysical objects in the MR plane, going from GCs to galaxies of different mass and morphological type and eventually to CoGs, suggests that a unique relation could exist for all of them and that such a relation likely owes its origin to the cosmological growth of DM halos. The distribution of the DM halos and their number density as a function of redshift has been the subject of several studies which culminated with the large-scale numerical simulations of the Universe. We cite here one for all, the Millennium Simulation (Springel et al., 2005). In parallel the studies of the halo growth function, HGF, as the integral of the halo mass function, HMF, appeared in literature [see, for instance, Lukić et al. (2007), Angulo et al. (2012), Behroozi et al. (2013)]. The HGF gives the number density of halos of different mass per (Mpc/ h 3 emerging at each epoch by all creation/destruction events and consequently yields the halos that nowadays populate the MR plane and generate the observed galaxies. Chiosi et al. (2019) adopted the HGF of Lukić et al. (2007) who, using the ΛCDM cosmological model and the HMF of Warren et al. (2006), derived the number density of halos n(M _D , z) over ample intervals of halo masses and redshifts. Since the n(M _D , z) of Lukić et al. (2007) refers to a volume of 1 (Mpc/ h 3 , before being compared with the observational data, it must be scaled by a suitable factor to match the volume sampled by observations. Anyway, the following characteristics of the HGF are worth being noted: 1) for each halo mass (or mass interval) the number density is small at high redshift, increases toward the present, and reaches a maximum at a certain redshift. The peak is either followed by a descent (for low mass halos) or a plateau (for high mass halos). In other words, first the creation of halos outnumbers the destruction, whereas the opposite occurs in general for low mass halos after a certain redshift. 2) At any epoch high mass halos are much less numerous than the low mass ones. This implies the existence of a cut-off mass at the high mass side. 3) The HGF also implies that halos of different mass have a given probability of existence at any redshift [see for more details Chiosi et al. (2012, 2019)].

Assuming a certain number density of halos N _s derived from the observational data, Chiosi et al. (2019) set up the equation n (M _D , z) = N _s whose solution yields the mass of the halos M _D (z) as a function of the redshift and vice-versa the redshift for each halo mass. In practice for any value N _s one gets a function M _D (z). To each value of M _D along this function, with the aid of Eqs. 11–12), one can associate a value of M _s and R _e . The MRR of luminous galaxies is the result.

Notably for the N _s corresponding to 10⁻² halos per (Mpc / h 3 (roughly the volume surveyed by the SDSS [see Chiosi et al., 2019, for details]), the curve R _e (M _s ) falls at the edge of the observed distribution of ETGs in the MR plane. Higher N _s would shift the curve to larger halos, the opposite for lower N _s . One can therefore draw in the MR-plane the locus of the most massive M _D and associated M _s imposed by the halo HGF. The equation n (M _D , z) = N _s with N _s = 10⁻² or equivalently 10⁶ halos per 10⁸ Mpc³ rewritten to derive the halo mass M _D as a function of z is log M D = 0.0031546 z 3 − 0.006455 z 2 − 0.183 z + 13.287 . (13)

Starting from this, Chiosi et al. (2019) associate M _s and R _e to each M _D for any value of the redshift. The best fit of the resulting MR relation, limited to the mass interval of normal ETGs, 9.5 ≤   log  M _s ≤ 12.5 (M _s in solar units), is log R e = 0.048562 ( log M s ) 3 − 1.4329 ( log M s ) 2 + 14.544 ( log M s ) − 50.898 . (14)

Note that 1) the locus on the MR-plane predicted by N _s = 10⁻² halos per (Mpc / h 3 nearly coincides with the observational MRR; 2) the slope gradually changes from 0.5 to 1 going from low masses to high masses in agreement with the observational data [see van Dokkum et al. (2010b), and references therein]; 3) finally, Eq. 14 is ultimately linked to the top end of the halo masses (and their associated baryonic objects) that might exist at each redshift. Chiosi et al. (2019) named this locus the Cosmic Galaxy Shepherd.

The extrapolation of the Cosmic Galaxy Shepherd downward to GCs and upward to CoGs yields the relation log R e = 0.007584 [ log ( m ⋅ M s ) ] 3 − 0.1874 [ log ( m ⋅ M s ) ] 2 + 1.908 [ log ( m ⋅ M s ) ] − 9.027 (15) where R _e and M _s are in the usual units and m is the ratio m = M _D /M _s , for which a mean value of m = 25 is adopted ² . As already said this equation represents the cut-off mass of the HDF at different redshift, however translated into the R _e vs M _s . This gives a profound physical meaning to the line splitting the MR-plane in two regions, i.e. the region where galaxies are found, and that of avoidance, the so-called Zone of Exclusion (ZoE) found by Burstein et al. (1997).

Along the Cosmic Galaxy Shepherd, cut-off masses and redshift go in inverse order: low masses (and hence small radii) at high redshift and vice-versa. More precisely, halos and their luminous progeny that are born (collapse) at a certain redshift and are now located along the theoretical MRR of Eq. 11 associated to that redshift. Along each MRR only masses (both parent M _D and daughter M _s ) smaller than the cut-off mass are in place, each of these with a different occurrence probability. Clearly the low mass halos are always more common than the high mass ones. We will argue that in the MR-plane, only the most massive GCs, DGs, and ETGs are expected to fall along the Cosmic Galaxy Shepherd. All other objects of lower mass, the DGs in particular, are expected to lie above this limit. This suggests that there are other physical processes concurring to shape the observed MRR. In other words, the question is “what really determines the position of each galaxy on the MR-plane?”

To answer the above question Chiosi et al. (2012, 2019) argue what follows. The gravitational collapse of a proto-cloud generating a luminous galaxy is surely accompanied by star formation, energy feed-back, gas cooling and heating, loss of mass and energy by winds, acquisition of mass and energy by mergers, etc. Therefore, the result of all these processes taking place together may largely differ from one case to another and also differ from the ideal case of a dissipation-less collapse. For this latter (Gott and Rees, 1975; Faber et al., 1984; Burstein et al., 1997) derived the relation R D ∝ M D 0.53 . (16)

Inside this halo a galaxy with stellar mass M _s and a half-mass radius R _e is built up over the years. Chiosi et al. (2019) take the dissipation-less collapse as the reference case. Using the data of the Illustris models, they derive the following MRRs log R e = 0.541 ⁡ log M s − 4.702 + k m f o r ⁡ log M s > 10.5 (17) log R e = 0.102 ⁡ log M s − 0.017 + k d f o r ⁡ log M s < 10.5 (18) where the constants k _m and k _d can be determined by fixing the initial conditions of the collapsing proto-halo. The slope of Eq. 17 does not significantly differ from that of the dissipation-less collapse, Eq. 6 and that of the empirical MRR of ETGs, Eq. 4. Along each MRR of the theoretical manifold, the agreement between data and theoretical models seems to be possible only for the most massive galaxies. For smaller masses, the slope of the theoretical MRR, Eq. 18, is much flatter than the observational one (about a factor of two).

From the above considerations one could suggest that the Cosmic Galaxy Shepherd and Eq. 4 represent the locus in the MR-plane of galaxies formed by quasi dissipation-less collapses. In contrast, special conditions ought to hold for all other objects that deviate from this condition. The explanation is different for the monolithic and hierarchical scenarios:

a) In the monolithic view, in addition to star formation, galactic winds are the key ingredient to consider, in particular for low mass galaxies, because DGs show the largest deviation from the observed MRR, Eq. 4 or Eq. 16. The analysis of the problem made by Chiosi et al. (2019) shows that: (i) the stronger the galactic wind the larger is the final R _e . Galaxies depart from the locus represented by Eq. 4 and/or Eq. 16 at decreasing mass and increasing galactic wind, the low mass ones having the strongest effect; (ii) the efficiency of winds tends to decrease at increasing initial density. This means that the inflating effect of galactic winds in low mass galaxies of high initial density is low and the final radius of these galaxies will be close to the value predicted by Eq. 4 and/or Eq. 16. In conclusion the flatter slope of the theoretical MRR is likely produced by galactic winds.

When does the MRR develop in the course of time and evolutionary history of galaxies? In Figure 5 we show the R _e vs M _s distribution of the Illustris models at four cosmic epochs. At high redshifts, the distribution is clumpy and irregular. However, starting from z ∼ 1.5 and more clearly at z = 0, a tail-like feature develops on the side of large masses, say for masses ≳ 2 ⋅ 10¹¹ M _⊙. The best fit at redshift z = 0, using the relationship log  R _e = ϵ  log  M _s + η (masses and radii are in M _⊙ and kpc), yields the following values: for log  M _s > 11.3 ϵ = 0.651 and η = − 6.557, while for log  M _s < 11.3 ϵ = − 0.005 and η = 0.592. What are the causes of the cloud-like and tail-like distributions? Why a cloud dominates the low mass range? Why the tail is well visible only for the high masses at low redshifts? Which is the physical meaning of this distribution? To cast light on this Chiosi et al. (2019) examined the history of R _e and M _s for several individual galaxies. The main conclusion of their analysis is that mergers among objects of low and comparable mass can generate galaxies with larger masses and radii, but exceptions are possible in which either the mass or the radius or both decrease. In general the galaxies do not leave the cloud region. All this does not contradict the previous case of Buonomo (2000) because the monolithic counterparts to compare with are not available. The cloud region is instead roughly coincident with the distribution of DGs of different types (see the discussion by Chiosi and Carraro, 2002). At the same time mergers among galaxies with different masses and/or comparable masses can generate objects that shift outside the cloud producing the MR-sequence (actually they define it), the locus of which agrees with the observed distribution for ETGs [see e.g. Chiosi et al. (2019), and references therein]. The stellar content of massive ETGs suggests that star formation has ceased long ago so that strong energy feedbacks are absent and the systems are close to the virial equilibrium. This implies that important mergers do not longer occur. At variance DGs are still undergoing frequent mergers, active star formation episodes, and strong galactic winds. They cannot be therefore in this ideal condition of equilibrium and so they depart from the observed MRR. Nevertheless, there are some DGs that fall along the MRR of massive ETGs and therefore are likely in a similar dynamical and star-forming condition, e.g. ωCen and M32 [see Chiosi et al. (2019), for more details].

On consideration of these premises, Chiosi et al. (2019) argued that the observed distribution of ETGs, inactive DGs and GCs, represents the locus of objects that have reached the ideal situation of mechanical equilibrium and pure passive evolution. They cannot go beyond this limit. Their MRR is therefore in the boundary between the permitted and forbidden regions of the MR-plane.

Putting the many tessarae of the mosaic together, the conclusion is that the observational MRR is the intersection of the theoretical manifold of the MRRs (each curve being labeled by the collapse redshift from the past to the present) with the Cosmic Galaxy Shepherd, along which objects in mechanical equilibrium and passive evolutionary state are located. To prove this statement Chiosi et al. (2019) resorted to the method proposed long ago by Chiosi and Carraro [see Chiosi and Carraro (2002)], however updating it with recent theoretical and observational data. In the MR-plane of Figure 4 they draw two loci and a mass interval as a function of the initial density (redshift):

1) The first locus is the MRR traced by models of different mass but same initial density and formation redshift. Using all models to disposal (Chiosi and Carraro, 2002; Merlin and Chiosi, 2006; Merlin and Chiosi, 2007; Merlin et al., 2010; Chiosi et al., 2012; Merlin et al., 2012; Chiosi et al., 2019), this locus is described by the relation log R e = [ − 1.172 − 0.412 ( 1 + z f ) ] + [ 0.244 + 0.0145 ( 1 + z f ) ] log M s . (19)

This expression is robust thanks to the regular behavior of the models and the density-mass-radius relationship of Eq. 16. Relation (19) is compatible with the MRRs predicted by Fan et al. (2010) and the models of Illustris by Vogelsberger et al. (2014). The cases shown in Figure 4 are: z _f ≃ 1, z _f ≃ 2, z ≃ 5 and z _f ≃ 10.

2) The second locus is the Cosmic Galaxy Shepherd. Among the various HGFs in literature (Lukić et al., 2007; Angulo et al., 2012; Behroozi et al., 2013), we adopt the HGF of Lukić et al. (2007) and make use of the analytical expression for the Cosmic Galaxy Shepherd extending across the whole MR-plane given by Eq. (15). However, to better illustrate this issue, we present here an analytical approach based on the classical halo mass distribution of Press and Schechter (1974) that is supposed to trace also the mass distribution of luminous galaxies (assuming one galaxy per halo). At each redshift, the HGF of Press and Schechter (1974) provides the relative number of galaxies per mass bin. The cut-off mass M D C O of the Press and Schechter (1974) function yields the maximum limit for the galaxy masses at each redshift. In the Press and Schechter (Press and Schechter, 1974) formalism, the cut-off mass varies with redshift according to: M D C O = M N × ( 1 + z ) − 6 n + 3 (20)

The exponent n represents the slope of the power spectrum perturbations and M _N is a suitable mass scale normalization. At any redshift, most galaxies have total masses smaller than M D C O , even if higher values cannot be excluded. It can be easily shown that the fractional mass in (or the fractional number of) galaxies with mass greater than M D C O is a function of n. For n = −1.8, the percentage of galaxies in the interval M D C O < M D < 10 × M D C O is about 15% while in the range 10 × M D C O < M D < 100 × M D C O is about 1%. Therefore, at any redshift galaxy masses up to say 50 × M D C O have a significant occurrence probability. Their radius is derived with the aid of the M _s vs M _D and R _e vs M _D relationships. For M D = γ M D C O , γ = 50, and n = − 1.8, one gets: R e = 16.9 × 1 0 − 12 × γ − 0.79 × ( 1 + z ) 3.96 × M s , (21) where R _e and M _s are in kpc and M _⊚. These are shown in Figure 4 with the dotted lines labeled by the redshifts z ≃ 1, ≃ 2, ≃ 5 and ≃ 10. On the MR-plane, they give the rightmost extension of the lines of constant density and hence they identify the maximum galaxy mass. At decreasing redshift this boundary moves progressively toward higher masses. Similar results can be obtained by means of the HGFs of Lukić et al. (2007), Angulo et al. (2012), Behroozi et al. (2013), the first of which is the Cosmic Galaxy Shepherd.

3) Finally, the third locus gives the expected interval for M _s for objects with total mass M _T between M D C O and 10 × M D C O as a function of redshift. Here the relation M _s (M _T ) has been plugged into Eq. 20 for M T C O . The permitted intervals are visible in Figure 4 by the horizontal lines labeled M s C O . The interval for M _s going from 10¹⁰ M _⊙ to 10¹² M _⊙ is fully compatible with the redshift interval for the formation of the majority of stars in a galaxy, i.e. from 2 to 1. This is also the mass range over which at any epoch the probability for the occurrence of massive galaxies falls to a negligible value. In different words, the right-hand border of the MRR has a natural width.

In this context, the relationship for ETGs, see Eq. 4, extended to the whole mass range from GCs to CoGs should correspond to the intersection between the lines of constant initial density and the lines where γ M T = M T C O ( z ) for equal values of the redshift (at least for all values of redshift > 1). This is what we see in Figure 4, i.e. the straight line marked by the large empty squares. This line nearly coincides with the Cosmic Galaxy Shepherd derived from the HGF of Lukić et al. (2007) that is marked by the crossed dark-red line in Figure 4, i.e. Eq. 14 and/or (15). Finally, this line is also coincident with the locus traced by objects that underwent a dissipation-less collapse (or very close to it) and are nowadays in mechanical equilibrium and passive evolutionary state. This is mainly traced by GCs, few DGs (the large majority of DGs lie above it), ETGs, and a number of CoGs. This confirms the result by D’Onofrio et al. (2020): only passive galaxies (strongly decreasing today in their luminosity) trace the MRR with a slope varying from 0.5 to 1, the highest value being reached by galaxies that suffered the strongest luminosity decrease with the redshift, i.e. those that long ago ceased their stellar activity, i.e. the most massive ones. Spirals occupy approximately the same location of ETGs in the MR-plane, thus suggesting that their ongoing star formation is not affecting the overall situation of mechanical equilibrium. Furthermore, it is worth noting that the slope of MRR derived from the HGF is about 1 in the range of massive galaxies (say above 10¹² M _⊙), i.e. formally identical to the MRR derived from the virial theorem. This coincidence might suggest a dependence of the observed MRR slope from the virial condition. The true driver is instead the HGF, more precisely its fall off toward high values of the halos’ masses at any redshift. To conclude, all the objects along the MRR are in virial conditions and passive evolutionary state (all mechanical process and star formation activity are at rest).

In the case of ETGs the best tool highlighting the main driver of the SF and the SFH is the NB-TSPH hydrodynamic simulations, in which the rate of star formation is usually expressed by the Schmidt (Schmidt, 1959b) law d ρ s d t = − d ρ g d t = c ∗ ρ g k t g (22) where ρ _s is the current mass density of stars, ρ _g is the current mass density of gas, t _g is a characteristic time scale (typically the free-fall), k is a suitable exponent (typically k ≃ 1), and c ^∗ is the so-called dimensionless efficiency of star formation (typically c ^∗ = 0.01 ÷ 0.1).

Based on simple arguments, there are at least three prerequisites for gas (likely in the form of molecular clouds) to be eligible to star formation: the gas has to be in convergent motion, i.e. the velocity divergence must be negative; the gas must be gravitationally unstable, i.e. it must satisfy the Jeans condition t _sound ≥ t _ff (where t _sound is the time scale related to the local sound velocity); the gas must be cooling, i.e. it has to verify the relation t _cool ≪ t _ff . Normally, SPH codes treat star formation simply implementing the Schmidt law in the computational language and transforming part of the gaseous particles that satisfy the three conditions above in new, collisionless particles of different mass (“stars”). The characteristic time scale is chosen to be the maximum between t _cool and t _ff time-scales; however in most situations t _g = t _ff is also a good choice. Knowing ρ _g and ρ _s and integrating upon the current volume of the system one gets the current values of M _g and M _s . Nowadays there are numerous galaxy models whose stellar content has been calculated with the above prescription. However they differ in a number of important assumptions, chief among others the cosmological model of the Universe and the scenario in which galaxy formation and evolution is framed. In the following, for the sake of illustration we will summarize here the results of three paradigmatic cases, i.e. the pure monolithic scheme of Chiosi and Carraro (Chiosi and Carraro, 2002), the early hierarchical scheme of Merlin et al. (Merlin et al., 2010; Merlin et al., 2012), and the full hierarchical scheme, e.g. the Illustris case of Vogelsberger et al. [e.g. Vogelsberger et al., 2014, and references]. It is worth recalling here that care must be paid on the link between the Schmidt and Kennicutt-Schmidt SF laws and their implications for numerical simulations (Schaye and Dalla Vecchia, 2008).

The pure monolithic scheme. Chiosi and Carraro (2002) highlighted the role of over-density of the initial perturbation when exceeding the threshold value. Two groups of models were analyzed according to the initial over-density: 1) models with mean initial density ⟨ρ⟩≃ 200 × ρ _u (z) and collapse redshift z _f = 5 (shortly named A); 2) models with ⟨ρ⟩≃ 5 × ρ _u (z) and z _f ≃ 1 (named B). ρ _u (z) is the density of the Universe at redshift z ³ . Perturbations with spherical symmetry, assigned mass, mean density exceeding the critical value (and hence suitable radius) are let collapse and form stars. A MonteCarlo procedure is adopted to fix the initial coordinates and velocities of the DM and BM particles. The key result of this study is that the star formation history (rate vs. time) is found to depend on the depth of the gravitational potential well of a galaxy. The following picture can be drawn. In the case of deep gravitational potentials (such as in massive and/or dense galaxies) once star formation has started energy is injected into the gas by supernova explosions, stellar wind, etc., but this is not enough to push the gas out of the potential well. The balance between cooling and heating is reached and the gas consumption by star formation goes to completion. Star formation cannot stop until the remaining gas is so little that any further energy injection will eventually heat it up to such high energies (temperatures) that the gravitational potential is overwhelmed. No more gas is left over and star formation is quenched. The star formation history resembles a strong unique burst of activity, a sort of monolithic star forming event, taking place over a certain amount of time, of the order of 1–2 Gyr. In contrast, in a galaxy of low mass and/or density and hence shallow gravitational potential, even a small star-forming activity will heat up the gas above the potential well. Some of it is soon lost in galactic wind, the remaining one becomes so hot that it will take long time to cool down and to form new stars. The cycle goes on many times in a sort of repeated bursting mode of star formation taking place during long periods of times if not forever. Out of all this we can derive what follows:

1) The duration, strength, and shape of the SFR as a function of time strongly depend on the galaxy mass and the initial density: (a) Galaxies of high initial density and total mass undergo a prominent initial episode of SF followed by quiescence. (b) The same happens to high mass galaxies of low initial density, whereas the low mass galaxies experience a series of burst-like episodes up to the present. The details of their SFH are very sensitive to the value of the initial density. The typical dependence of the SFR on time for models B is shown in the left panel of Figure 6, while the right panel shows the SFR of a low mass galaxy (M _T = 10⁹ M _⊙) for moderate variations of the initial density. Models A are not shown because their SFR is represented by a single initial spike.

The early hierarchical scheme. Merlin et al. (2010, 2012) using initial conditions derived from large-scale cosmological simulations and abandoning the strict monolithic scheme, much improved the NB-TSPH galaxy models of Chiosi and Carraro (2002). They also adopted the ΛCDM cosmology instead of the classical CDM ⁴ . Without entering into detail, they cut from large-scale simulations calculated with the free code COSMICS by Bertschinger (1995, 1998) a spherical portion containing proto-halos of DM and BM in cosmological proportions with the desired over-density, mass, and size (the reference proto-halo with the highest mass to consider). The same is made for halos with lower mass and smaller dimensions at fixed mean density. The procedure to obtain halos with the same mass but different initial mean density is more complicated and will not be reported here [see Merlin et al. (2010, 2012) for the details]. These proto-halos contain a number of distinct lumps of matter that will merger together later on. The cosmological simulation provides the initial positions and velocities of all the particles in the proto-halos. The expansion of the Universe is taken into account. The proto-halos are followed through their expansion (caused by the Hubble flow), down to their collapse and aggregation into single objects. The redshift of the collapse varies from model to model and, inside the same model, from the center to the periphery. In general the collapse occurs in the redshift interval 4 > z > 2, it starts in the central regions and gradually moves outward. The collapse is complete at redshift z ≃ 2. All models develop a stellar component. The more massive halos experience a single, intense burst of star formation (with rates ≥ 10³ M _⊚/yr) at the early epochs. The intermediate mass halos (M _T ≃ 10¹¹ M _⊙) have star formation histories that strongly depend on the initial over-density, i.e. with a single or a long lasting period of activity and strong fluctuations in the rate. The small mass halos (M _T ≃ 10⁹ M _⊙) always have fragmented star formation histories: this is the so-called galactic breathing phenomenon. These models are classified as early hierarchical because they experience repeated episodes of mass accretion at very early epochs and then evolve in isolation ever since. They confirmed the correlation between the initial properties of proto-halos and the star formation history found by Chiosi and Carraro (2002). The models have morphologies, structures, and photometric properties similar to real galaxies [see Merlin et al. (2012),  for all other details].

The fully hierarchical scheme. This is the most difficult case to discuss because of mergers among galaxies of different mass, size, and age. If gas is present recurrent episodes of stronger star formation activity may occur. It is conceivable that seed galaxies prior to any encounter behave like the general scheme envisaged before and governed by the initial density and mass. Mergers among objects of similar mass would likely enhance the rate of star formation in a sort of burst of short duration and fold the two histories together. Mergers among objects of much different mass would simply generate a temporary perturbation on the star formation history of the most massive one, while less massive object simply loses its identity. In the hierarchical scenario, tracing the star formation history of single galaxies is a hard task. Anyway, the observational evidence provided by the stellar content of galaxies of different mass strongly supports the mass-density scheme we have described.

The general trends of the SFR described in this section agree with the picture envisaged long ago by Sandage (1986), examining the SFR in galaxies of different types [see also Tinsley (1980), Chiosi et al. (2014), Matteucci (2016)]. This scenario has been confirmed by studies of SF histories based on absorption line indices (Thomas et al., 2005), and by the recent study of Cassarà et al. (2016). A good agreement also exists with other independent numerical NB-TSPH models of galaxy formation and evolution by Kawata and Gibson (2003a, 2003b), Kobayashi (2005).

According to Matteucci (2016) the most common parameterization of the SFR in LTGs is the Kennicutt (Kennicutt and Jr., 1998a) generalization of the original Schmidt (Schmidt, 1959b) law, where the SFR is proportional to the gas density ρ. Kennicutt (Kennicutt and Jr., 1998a) suggested that the SFR can be written as: S F R ( t ) = ν Σ g κ (23) where Σ is the gas surface mass density, ν the efficiency of star formation (the SFR per unit mass of gas), and κ = 1.4 ± 0.15, as deduced by the data of the star-forming galaxies [see also Kennicutt (1998)]. Other parameters, such as gas temperature, viscosity, and magnetic field, are not considered.

Actually, the “Kennicutt law” was in use long before its discovery. In the mid-seventies Larson (1975, 1976) developed the first modern hydrodynamic models of formation and structure of elliptical and spiral galaxies, showing that a rate of star formation strongly declining during the latest stages of collapse was necessary to form a massive disk in spiral galaxies. However, once the gas has settled onto the equatorial plane and built up the disk, the rate of star formation should increase to a peak value and then decline again. The duration of this phase and the height of the peak were found to depend on the position on the disk. Larson envisaged several physical mechanisms that might strongly suppress star formation during the latest stage of collapse, e.g. velocity dispersion of the gas, tidal forces exerted on the remaining gas by the already formed spheroidal component, and dependence on the cloud-cloud collision frequency. The same processes were also invoked to control the second phase of star formation. Starting from this Talbot and Arnett (1975) correlated the process of star formation with the surface mass density of the gas in an already flattened disk, whose thickness is regulated by the balance between the gravitational attraction and the increase of the scale height by energy injection by short-lived stars (e.g. type II supernova explosions by massive stars). They proposed a star formation rate proportional to the surface mass density of gas. Chiosi (1980) folded the Larson (Larson, 1976) results into the Talbot and Arnett (Talbot and Arnett, 1975) mechanism and incorporated all this into a new model for the chemical evolution of galactic disks in the presence of infall. In this model the disk is described by a series of concentric rings (no mass exchange among them), whose surface mass distribution at the present time t _g is given by an exponential law of type Σ(r) = Σ _d exp( − r/R _d ), where Σ _d and R _d are two scale parameters. The formation of the disk is supposed to occur by rapid infall of the gas left over by the formation of the halo and the central spheroidal component. The temporal and spatial dependence of the infall rate is given by d Σ ( r , t ) d t = A ( r ) e x p ( − t / τ ) (24) where A(r) is a suitable function to be determined. This is derived by integrating Eq. 24 with respect to time and by equating it to the present-day mass distribution. We obtain A ( r ) = Σ d e x p ( − r / R d ) τ − 1 [ 1 − e x p ( − t g / τ ) ] − 1 (25) for r _B ≤ r ≤ r _D , where R _B ≃ 2 kpc and r _D ≃ 20 kpc are the typical radius of a bulge and of a disk of spiral galaxies, respectively. The scale parameters Σ _d and R _d are determined by knowing the rate and the surface mass at a certain position of the disk (e.g. the solar vicinity in our case). Thanks to the short time scale of the energy input from massive stars (a few million years), compared to the mass accretion time scale by infall (from hundred to thousand million years) the disk was supposed not to differ from an equilibrium state so that the Talbot and Arnett (Talbot and Arnett, 1975) formalism could be applied. Chiosi (1980) and Chiosi and Matteucci (1980) proposed and used the SFR: d Σ s ( r , t ) d t = − d Σ g ( r , t ) d t = ν ̃ Σ ( r , t ) Σ g ( r , t ) Σ ( r ̃ , t ) κ − 1 Σ g ( r , t ) (26) where Σ _g (r, t) and Σ _s (r, t) are the surface mass densities of gas, stars at the position r or and time t, respectively. The quantities Σ ( r ̃ , t ) and ν ̃ are the total surface mass density at a particular distance from the galaxy center, and an efficiency parameter. They play the role of a particular radial scale controlling star formation. In the Larson’s view they might be associated with the radial distance at which the central spheroidal component and the innermost regions of the disk exert their tidal effect on the residual external gas. The spatial and temporal dependence of the relation Eq. 26 in the infall model for the disk of the Milky Way and disk galaxies in general is such that at any time the SFR is strongly inhibited at distances r > r ̃ , while at any given r the SF starts small, increases to a peak value, and then declines again. This behavior of the SFR is typical of all infall models, where because of interplay between gas accretion and consumption, the SFR starts low, reaches a peak after a time approximately equal to τ, and then declines. Independently of the position, the net temporal dependence of the SFR is the time delayed exponentially declining law: S F R ∝ t τ exp − t τ . (27)

The Schmidt law is the link between gas accretion by infall and gas consumption by star formation. Thanks to the infall model by varying τ (time scale of the galaxy formation process) one can recover all types of star formation indicated by observational data going from GCs to LTGs and ETGs. The infall scheme and companion SFR have been widely used in many studies on the subject of galactic chemical evolution [e.g. Matteucci (2016), for a recent review and references]. The infall galaxy model is very flexible and can be adapted to a wide range of astrophysical problems. Suffice to recall that it has been used by Bressan et al. (1994) to model the spectro-photometric evolution of ETGs reduced to point mass objects, extended by Tantalo et al. (1998a) to the case of spherical systems made of BM and DM mimicking ETGs, adapted by Portinari and Chiosi (2000) to include radial flows of gas in disk galaxies, and recently used by Chiosi et al. (2017) to study the cosmic star formation rate and by Sciarratta et al. (2019) to investigate the color-magnitude diagram of galaxies in general.

The connection between the structure and dynamics of galaxies and their stellar population, which we have encountered addressing the FP problem, is also part of the SF problem of galaxies. Observations have in fact revealed that the SFR and the stellar mass (M _s ) of active star-forming galaxies are tightly correlated ( S F R ∝ M s 0.6 ). This trend is known as the galaxy “main sequence” (MS) (Brinchmann et al., 2004; Noeske et al., 2007; Salim et al., 2007; Elbaz et al., 2011; Salim et al., 2012; Whitaker et al., 2012; Rodighiero et al., 2014; Speagle et al., 2014; Schreiber et al., 2015). Adopting different samples the MS may be different, either in slope or in scatter, primarily for selection effects on the adopted SF indicator used (Popesso et al., 2019). One may select galaxies according to their mass and/or color, picking preferentially the blue cloud objects, or using the BzK color selection (Daddi et al., 2004; Daddi et al., 2007; Pannella et al., 2009) or the UVJ selection (Williams et al., 2009; Whitaker et al., 2012), or adopting a minimum threshold for the specific SFR (sSFR = SFR/M _s ) (Karim et al., 2011).

The presence of a main sequence, with a scatter of 0.3 (in log units for active star-forming objects), indicates that these galaxies have an SFR that spans a factor of two. This can be explained by the self-regulating nature of the SF process, that is by the interplay between gas accretion, SF, and feedback (Schaye et al., 2010; Davé et al., 2011; Haas et al., 2013; Lilly et al., 2013; Rodríguez-Puebla et al., 2016; Tacchella et al., 2016).

The scatter however is much larger ( ∼ 0.6) if all types of galaxies are considered. We can see it in Figure 7. The red dots in the various panels represent the data of the WINGS database (Fritz et al., 2007; Fritz et al., 2011). The SFR in the last 20 Myrs, measured from the spectral energy distribution in more than 3,000 objects of all morphological types, is plotted versus the stellar mass. The artificial data coming from the Illustris simulations are also represented for different redshift epochs: z = 4 (blue dots), z = 1 (green dots), and z = 0 (black dots). Before drawing any conclusion, it is worth recalling that the model galaxies of the Illustris simulation were chosen to have stellar masses above 10⁹ M _⊙ at z = 0. Therefore, the comparison between theory (black dots) and data (red dots) in Figure 7 is possible only for M _s ≥ 10⁹ M _⊙. We note that in the common region the simulations predict the correct slope and quite a similar scatter at z = 0. They also predict that the slope mildly changes with redshift and that the scatter increases going to the present epoch. This limit on M _s does not exist for the samples at higher redshifts.

The origin of the scatter and its amount might be different for dwarf and giant galaxies (Matthee and Schaye, 2019) and can be attributed to both short-time and long-time processes, such as the competing effects of inflows and outflows, the variation of the halo mass, the variation of the SFE, and the feedback effects from the active nuclei. Daddi et al. (2007), Elbaz et al. (2007), Noeske et al. (2007) claim that the correlation is present even at redshift ∼ 2, with a nearly constant slope and a dispersion similar to that observed for galaxies in the local Universe (Brinchmann et al., 2004).

The information that one can draw from the MS is still under debate (Kelson, 2014; Abramson et al., 2015). There are many open questions: does the main sequence imply that the SFH of galaxies of the same stellar mass is similar? Is the MS a median “attractor-solution” (Peng et al., 2010; Behroozi et al., 2013)? Is it an average sequence for a population at a certain age of the Universe (Gladders et al., 2013; Abramson et al., 2016)? Do galaxies of the same mass have different SFHs on longer time-scales? What effects are most significant at different mass and time-scales?

The slope and scatter of the MS might encode such crucial information. What makes the growth of galaxies different? Which are the important time-scales of SF? Which are the systematic and stochastic effects behind the scatter? Which is the role played by the environment and by the DM in the assembly accretion history?

The MS measured at higher redshifts shows a positive correlation evolving only a bit in slope and scatter [see e.g. Elbaz et al., 2007, Noeske et al., 2007, Whitaker et al., 2012]. This might support the idea that the link between structure and stellar population in galaxies is already in place at z ∼ 2.5 (Wuyts et al., 2011).

The SFR in a galaxy depends on a variety of factors, such as the rate at which the galaxy accretes mass from the IGM, the rate of shocking and cooling of gas onto the galaxy, the details of how the inter-stellar medium (ISM) converts gas into stars, the amount of galactic fountain and outflow, etc. This complex nonlinear physical mechanism is difficult to understand, in particular if one wants to discover what processes dominate, and if and how these change over time. The PCA reveals that neutral gas fraction f _gas , stellar mass M _s , and SFR form a nearly flat 2D surface (Lagos et al., 2016). The location of the plane varies with redshift, and galaxies can move along it when f _gas and SFR drop with redshift. Their position along the plane is correlated with gas metallicity. This is a sort of “fundamental plane of SF” whose curvature is determined by the dependence of the SFR on gas density and metallicity.

The previous sections have clearly demonstrated that the observed properties of galaxies are regulated by a complex series of physical effects tightly intertwined. Last, but not the least, is the ratio between the stellar mass in a galaxy and its dark matter component M _s /M _D (and its inverse M _D /M _s ). The ratio M _s /M _D is a quantity that ultimately affects the half-luminosity radius R _e of the stellar component of a galaxy (see Section 6.1). The analysis of the Illustris data and the theoretical galaxy models of Chiosi and Carraro (2002), Merlin and Chiosi (2006), Merlin and Chiosi (2007), Merlin et al. (2010), Merlin et al. (2012), Chiosi et al. (2012) led Chiosi et al. (2019) to suggest that the ratio M _s /M _D depends on the total mass of the galaxy M _T ≃ M _D and the redshift z _f at which the bulk of SF occurs. This is shown by Figure 8 for the Illustris data. For low values of the redshift (say below 0.6), the ratio smoothly decreases with mass M _D (low mass galaxies are slightly more efficient in building their stellar content); the opposite occurs for higher redshifts, where M _s /M _D increases with M _D . Chiosi et al. (2019) give the following analytical expression for the ratio M _s /M _D as a function of M _D and z log M s M D = [ 0.218 z − 0.101 ] log M D + [ 0.169 z − 2.227 ] (28) where the halo mass goes from 10⁴ M _⊙ to 10¹⁴ M _⊙ and the redshift from 0 to 4. The ratios M _s /M _D predicted by Eq. 28 are indicated by the small black dots of Figure 8.

Other relationships for the inverse ratio m = M _D /M _s can be found in the literature [see, for instance, Shankar et al. (2006), Fan et al. (2010), Girelli et al. (2020)]. For M _D ≥ 10¹¹ M _⊚ Fan et al. (2010) propose the relation: m = M D M s = 25 M D 1 0 12 0.1 1 + z 4 − 0.25 (29) from which one derives the ratio M _s /M _D shown in Figure 8 by the red circles. In practice there is no dependence on redshift.

Notably, the curve of Fan et al. (2010) agrees with the one derived by Chiosi et al. (2019) using the Illustris models for the values of the redshift smaller than about 1.6 (the slope is nearly identical). Shankar et al. [Shankar et al. (2006), and references therein] presented a detailed analysis of the dependence of M _s on M _D . First, they claim that for M _D < 10¹¹ M _⊙ the relation should be m = M D M s = C M D − 2 / 3 (30)

with C a suitable proportionality constant to be determined. Assuming equality between the values of m derived with the two above relationships (at the transition mass M _D ≥ 10¹¹ M _⊚), the proportionality constant is log  C = 9.044. The ratios M _D /M _s resulting by Eq. 30 are shown in Figure 8 with the dark golden circles. Note that the relation of Shankar et al. (2006) agrees with that of the Illustris models for redshifts in the range from 2 to 4.

It is also worth noting that the linear extrapolation of the Fan et al. (Fan et al., 2010) relationship (red circles) in Figure 8 to lower masses and the linear extrapolation of the Shankar et al. (Shankar et al., 2006) curve (dark golden circles) to higher values of the mass encompass the predictions derived from the Illustris models for all the values of the redshift.

Shankar et al. (2006) derived a second analytical expression for the relation between M _s and M _D : M s = 2.3 × 1 0 10 M ⊚ ( M D / 3 × 1 0 11 M ⊚ ) 3.1 1 + ( M D / 3 × 1 0 11 M ⊚ ) 2.2 (31) for M _D ≥ 10¹¹ M _⊙. In this relation there is not an explicit dependence on the redshift. The ratios M _s /M _D predicted by Eq. 31 are visible in Figure 8 with the black filled squares. Eq. 31 predicts ratios m(M _D , z) that are in agreement with those from Eq. 28 derived from the Illustris data, Eq. 29 from Fan et al. (2010), and Eq. 30 only in the region around log(M _D ) ≃ 12 and z ≃ 0.

In a very recent study Girelli et al. (2020) have thoroughly investigated the stellar-to-halo mass ratio of galaxies (M _s /M _D ) in the mass interval 10¹¹ < M _D < 10¹⁵ and redshifts from z = 0 to z = 4. They use a statistical approach to link the observed galaxy stellar mass function on the COSMOS field to the halo mass function from the ΛCDM-Dustgrain simulation and derive an empirical model to describe the variation of the stellar-to-halo mass ratio as a function of the redshift. Finally they provide analytical expressions for the function M _s (M _D , z). The relations M _s /M _D vs M _D as a function of the redshift obtained with the formalism of Girelli et al. (2020) are also shown in Figure 8 (the magenta and dark-olive-green open circles joined by dashed lines of the same color). See also for a similar analysis the study of Engler et al. (2020).

It is soon evident that while all studies agree on the M _s /M _D ratios for objects with halo mass in the interval 11.5 ≤ logM _D ≤ 12.5 nearly independently of the redshift, they badly disagree each other going to lower values of the halo mass. Furthermore, they also disagree with the theoretical results predicted by Illustris. The problem is open to future investigations.

When galaxy formation started DM and BM were in cosmological proportions (i.e. M _D = ωM _B with ω ≃ 6). Then the SF gradually stored more and more BM into stars.

Here, exploiting again the Illustris library of model galaxies (Vogelsberger et al., 2014) we show the relationships between the stellar mass M _s (as a proxy of the BM component) and the dark mass M _D , and that between R _e and R _D for four different values of the redshift (z = 4, 2, 1, and 0). They are visible in the left and right panels of Figure 9, respectively. Masses (in M _⊚) and radii (kpc) are in log units and the color code indicates the redshift (z = 4, blue; z = 2, green; z = 1, yellow; z = 0, red).

It is clear that the efficiency of SF over the Hubble time, i.e. the transformation of gas in stars, is different in galaxies of different masses. Since M _B < M _D , M _s is always smaller than M _D . However, galaxies of different total mass can build stars at different efficiencies, and the ratio M _s /M _D is therefore expected to vary with M _D and redshift. In the left panel of Figure 9, we note that M _s increases with M _D , so that low mass galaxies build up less stars than the more massive ones. The slope of the relation however decreases as the redshift goes to zero. In more detail, for redshifts z ≳ 2 and masses M _D ≃ 10¹² M _⊚ the slope decreases at decreasing redshift so that more and more stars are present at given M _D . More precisely, for z ≲ 2 and M _D ≤ 10¹² M _⊙ the above trend holds, but above this limit the opposite occurs, at a given M _D less stellar mass is present than expected. In other words, massive galaxies are less efficient builders of their stellar content. We can approximate this relation between   log  M _s and   log  M _D with the linear dependence   log  M _s = α  log  M _D + β, where α and β may vary with the mass range and the redshift. From the linear fit we obtain: (z = 4, α = 1.55, β = − 8.19) (z = 2, α = 1.44, β = − 6.78) (z = 1, α = 1.16, β = − 3.37, for M _D < 12.0) (z = 1, α = 0.76, β = − 2.30, for M _D > 12.0) (z = 0, α = 0.93, β = − 0.43, for M _D < 11.5), and (z = 0, α = 0.79, β = − 1.22, for M _D > 11.5). The ratio M _s /M _D varies from 0.2 to 0.05 when the mass M _D increases from 10⁷ to 10¹² M _⊚ with mean value ≃ 0.10. The overall process of star formation is not highly efficient, large amounts of gas remain unused and likely expelled into the external medium through galactic winds partially enriched in metals. Similar results are given by Merlin et al. [Merlin et al. (2012), and references]. The efficiency of star formation is customarily measured by the ratio M _D /M _s as a function of M _D . This is simply given by: M D M s = 1 0 − β M D 1 − α (32) that has already been discussed in Section 11.

Similarly we can derive the relations:   log  R _e = γ  log  R _D + η ( R e = η R D γ ) that are shown in the right panel of Figure 9. From the linear fit we obtain: (z = 4, γ = 0.39, η = − 8.19) (z = 2, γ = 0.30, η = − 6.78) (z = 1, γ = 0.22, η = − 3.37, for M _D < 12.0) (z = 1, γ = 0.22, η = − 2.30, for M _D > 12.0) (z = 0, γ = 0.29, η = − 0.43, for M _D < 11.5), and (z = 0, γ = 0.29, η = − 1.22, for M _D > 11.0).

The radius R _D is larger than R _e by a factor of 3–10 as the galaxy mass increases from 10⁹ M _⊙ to 10¹³ M _⊙. The slope γ of the R _e -R _D relation (in log units) first decreases by about a factor of 2, from z = 4 to z = 1, and then increases again at z = 0. What is important is that while at high redshifts (our z = 4, z = 2, and z = 1 cases) the galaxy distribution on the R _e -R _D plane is a random cloud of points, at z = 0 a regular trend appears and R _e increases with R _D on the side of large values of R _D (largest masses). However, in the region of low radii and masses a cloud of points is still visible. The reason must be attributed to the effect of strong galactic winds and mergers among galaxies of similar mass in the hierarchical process that strongly perturb the mechanical equilibrium of these systems [see Chiosi et al. (2019)]. Finally, note that the ratios M _s /M _D ≃ 0.1 and R _e /R _D ≃ 0.1 − 0.3 confirm the predictions of Bertin et al. (1992), Saglia et al. (1992) based on analytical models for galaxies made of DM and BM.

The method employed for modeling stellar system in dynamical equilibrium is that of orbit superposition (Schwarzschild, 1979). The gravitational potential is defined as the sum of the central black hole (assumed a central point whose mass is to be determined) and of the stellar mass density derived from the stellar mass-to-light ratio. What is computed is the combination of orbits compatible with the spatially resolved stellar kinematics and photometric profiles. For the kinematically hot galaxies the early way to get the BH mass was based on the fit of the line-of-sight velocity dispersion of spherical galaxies assuming that the stellar distribution function is isotropic (Young et al., 1978). In more modern approaches, the fit is made over the entire line-of-sight velocity distribution (Rix et al., 1997; Gebhardt et al., 2007) for arbitrary galaxy models whose gravitational potential includes the effect of dark matter, and of triaxiality (Gebhardt and Thomas, 2009; van den Bosch et al., 2008). The most general and accurate possible models, with the highest resolution of spectroscopic observations, are reputed to be most accurate (Kormendy and Ho, 2013), provided that the BH sphere of influence is adequately resolved. A case in point is the estimate of the black hole mass in M87: early estimates yielded a mass ∼ 5 × 10⁹ M _⊚ from spherical, isotropic models Young et al. (1978). More recent analyses based on stellar dynamics yielded M _BH in the range ≈ (6. − 6.5) ⋅ 10⁹ M_⊚ (Gebhardt and Thomas, 2009; Gebhardt et al., 2011). The stellar dynamics mass value has been spectacularly confirmed by the Event Horizon Telescopes observations that yielded M _BH ≈ 6.5 ± 0.2|stat ± 0.7|sys ⋅ 10⁹ M_⊙ from the inference of the angular size of the black hole gravitational radius (Akiyama et al., 2019).

One of the most promising developments in the last years has been the increasing number of dynamical mass estimates obtained with ALMA [e.g., (Barth et al., 2016; Boizelle et al., 2019), for mildly active and quiescent galaxies]. ALMA has the capability to resolve cold molecular gas kinematics on angular scales well below 1 arcsec. (Wootten and Thompson, 2009). This is becoming instrumental to high-precision measurements of black hole masses in the “intermediate” mass domain, a previously uncharted territory. For instance, sub-parsec resolution ALMA observations revealed a black hole with mass ∼ 5 ⋅ 10⁵ M_⊚ in the dwarf galaxy NGC404 (Davis et al., 2020).

Space-based long-slit spectra of optical emission lines yield a velocity cusp (Macchetto et al., 1997). A striking example is provided by the radial velocity curve of NGC4374 (Bower et al., 1998): the STIS spectra show a Keplerian swing beginning at ± 0.5 arcsec and culminating at ± 0.1 arcsec, with a radial velocity difference of δv _r ≈ 400 km s⁻¹, implying an M _BH ≈ ( 1 . 5 − 0.6 + 1.1 ) ⋅ 1 0 9 M_⊚. The main concern is that gas motions could be affected by radiation forces, shocks, turbulence, and magnetic fields, and not only by gravitation. Relatively few galaxies have been found to have regular disk-like profile suggestive of a velocity field dominated by Keplerian motion in a dynamically cold disk (Kormendy and Ho, 2013). In addition, the Keplerian assumption is not consistent with gas flow toward low-accretion-rate SMBHs and at variance with observations of the Galactic Center (Jeter et al., 2019). For M87, both an early and a more recent analysis based on HST data suggest a black hole mass of ≈ (3 − 3.5) ⋅ 10⁹ M_⊚ (Harms et al., 1994; Macchetto et al., 1997; Walsh et al., 2013), and very close to the value obtained by modeling the jet boundary shape (Nokhrina et al., 2019), but always at variance with the values obtained from stellar dynamics (Section 13.1).

As mentioned earlier, the correlation between M _BH and host galaxy bulge properties—M _bulge and σ _⋆ or even bulge luminosity—is now an established fact since more than 20 years [see e.g. Kormendy and Richstone, 1995, Ferrarese and Merritt, 2000, Gebhardt et al., 2000]. Widely used forms of the relation between M _BH and σ _⋆ based on sources for which there is a dynamical M _BH determination are the ones of McConnell et al. (2011) and of Kormendy and Ho (2013) for early-type galaxies. Kormendy and Ho (2013) derived a power law: log ( M BH ) ≈ 8.491 ± 0.049 + ( 4.384 ± 0.287 ) log σ ⋆ , 200 , (34) where the mass is in solar units and σ _⋆ is units of 200 km s⁻¹. For both early type and spiral galaxies McConnell and Ma (2013) yield a significantly steeper slope ≳ 5, with a lower intercept for spiral galaxies, implying that M _BH in ETGs is about a factor 2 higher than in LTGs at a given σ _⋆ (McConnell and Ma, 2013). Equivalent relations (i.e., with similar scatter, around 0.30 dex) have been defined with the bulge mass, and infrared luminosity, usually suggesting a power-law relation between M _BH and M _bulge with an exponent ≈1 or larger (Kormendy and Ho, 2013; McConnell and Ma, 2013; Bennert et al., 2021).

There is still much ongoing research considering the linearity of the relation, its slope, and the origin of its dispersion. Theories that connect galaxy evolution and black hole growth predict the existence of a second parameter which may account for the dispersion in the M _BH—σ _⋆ correlation. Black hole—spheroid coevolution models would require that the BH mass scales with the gravitational binding energy of the spheroid host, ∼ M _bulge/r (Hopkins et al., 2007). The correlation can be easily turned into bivariate relations M BH ∝ M BH 0.6 σ ⋆ 1.2 and ∝ r bulge 0.6 σ ⋆ 2.4 that imply correlations between M _BH and σ _⋆ and M _bulge consistent with the observed ones (Saglia et al., 2016). A correlation with the binding energy of the host galaxy (Aller and Richstone, 2007; Barway and Kembhavi, 2007) implies the presence of a second parameter that may compensate for the changes in the galaxy structural parameters occurring at fixed M _BH. Saglia et al. (2016) found significant bivariate correlations consistent with a connection between M _BH and binding energy and with bulge kinetic energy, although the scatter remains comparable to the one for the M _BH—σ _⋆ correlations obtained with the best dataset of M _BH dynamical mass measurements.

The log(M _BH ) −   log(M _s,sph ) relation reported by Kormendy and Richstone (1995), Franceschini et al. (1998), Magorrian et al. (1998), McLure and Dunlop (2002), Marconi and Hunt (2003), Hring and Rix (2004) is almost linear, but the inclusion of low-mass spheroids revealed departures from linearity. Laor (1998, 2001), Wandel (1999), and Ryan et al. (2007) obtained a much steeper power law with a slope of 1.53 ± 0.14. The mean M _BH /M _s,sph ratio is probably not a universal constant, as it drops from ∼ 0.5% in bright (M _V ∼ − 22) ellipticals to ∼ 0.05% in low-luminosity (M _V ∼ − 18) bulges. Salucci et al. (2000) claimed that the M _BH − M _s,sph relation is significantly steeper for spiral galaxies than for (massive) elliptical galaxies. Graham (2012) suggested that the relation between luminosity (L) and stellar velocity dispersion (σ) for low-luminous ETGs is inconsistent with the M _BH − L and M _BH − σ relations. They prefer a broken M _BH − M _s,sph power-law relation, with a near-linear slope at the high-masses and a near-quadratic slope at the low-masses. In a recent review article Graham (2016) analyzed the consequences of this steeper relation, which can be rich of implications for the theories of galaxy–BH coevolution. Scott et al. (2013), Graham and Scott (2013) offered an interpretation for the curvature of the M _BH − M _s,sph relation invoking the presence of core-Sérsic and Sérsic spheroids at the high- and low-mass ends of the distribution respectively.

The highest degree of correlation is obtained for ETGs and for bulges that follow a Sérsic—surface brightness profile. Galaxies obeying a Sérsic photometric profile down to the resolution limits of their surface brightness profiles are believed to be the product of wet mergers, i.e., merger of gas rich galaxies that provide material to sustain accretion on the central black hole and trigger a period of sustained nuclear activity. The ensuing feedback effects (both radiative and mechanical) on the host, due to the active nucleus and to the merger-induced star formation, make it possible to couple the growth of the central black hole to the host spheroid mass (Zubovas and King, 2012; King and Pounds, 2015): the feedback forces by the quasar expel so much gas to quench both star formation and stop black hole growth, ultimately accounting for the relation between M _BH and σ _⋆, and M _BH and M _bulge (Di Matteo et al., 2005; Robertson et al., 2006).

However, the most massive elliptical galaxies often exhibit surface brightness profiles that are flatter than the extrapolation of Sérsic-like profiles. Sources showing a deficit with respect to the Sérsic profile are contributing to the scatter in the M _BH—M _bulge relation (Kormendy and Bender, 2009). Core profiles are believed to be due to dissipationless mergers of galaxies that have central black holes. N-body simulations show that merging of two galaxies with a sharp cusp may result in a merger remnant with a shallower core (Milosavljević and Merritt, 2001; Kulkarni and Loeb, 2012; Bortolas et al., 2016). The formation of a core has been ultimately linked to a bound binary black hole system, which produces a depletion of the stellar component in the nucleus due to slingshot ejection of stars on nearly-radial orbits.

The size of the core and the starlight and mass deficits in the centers of core galaxies (i.e., the mass ejected by the binary) have been found to scale approximately with the mass of the central black hole (Graham, 2004; Kormendy and Ho, 2013), in agreement with theory that predicts a mass deficit (Merritt, 2006) to be 0.5 n M _BH, with n as the number of major merger events. The luminosity deficit correlation provides an independent way to estimate M _BH in core ellipticals. Core radius is most strongly correlated with the black hole mass and correlates better with total galaxy luminosity than it does with velocity dispersion (Rusli et al., 2013). In addition, core scouring changes the orbit distribution. Only radial orbits allow for close passage past the galaxy center and thus only those stars can reach the vicinity of the central binary black hole. Consequently, the orbital structure in the core after core scouring is predicted to be strongly biased in favor of tangential orbits, while the ejected stars contribute to enhanced radial motions outside the core (Quinlan and Hernquist, 1997; Milosavljević and Merritt, 2001). For example, the orbital structure of the S0 NGC524 shows tangential anisotropy right at the SMBH radius of influence, corresponding to the core region in the photometric profile (Krajnović et al., 2009). Similar results apply to the elliptical galaxy NGC1600 (Thomas et al., 2016), and agree well with predictions from numerical simulations where core profiles are the result of SMBH binaries impoverishing the central nuclear regions (Rantala et al., 2018).

Recent work emphasizes the presence of substructures in the M _BH—σ _⋆ relation (Sahu et al., 2020). Pseudo-bulges are associated with spiral galaxies, and studies of their photometric profiles reveal that they are disk-like with a different surface brightness profile than classical bulges (Kormendy and Kennicutt, 2004; Kormendy and Bender, 2012b). Pseudo-bulges are known to be offset in the M _BH—σ _⋆ relation in the sense of systematically lower M _BH (Saglia et al., 2016). In the case of pseudo-bulges, the growth of the central black hole may be decoupled from the growth of the host spheroid and not associated with galaxy merger, but instead with mechanisms of secular evolution not related to gravitational interaction with other galaxies; in observational terms, some studies (Kormendy and Ho, 2013) find weak M _BH correlations for pseudo-bulges, see, however, e.g., (Bennert et al., 2021).

The most massive BHs have been detected only in the more luminous galaxies ( − 22 ≤ M _B ≤ − 18) (Ferrarese and Ford, 2005) and it is not clear yet if fainter and less massive systems host massive BHs and whether they follow the extrapolations of the SRs defined by the brightest objects. Searches for BHs in less luminous galaxies of the Local Group have produced ambiguous results, as in the case of M33 (Gebhardt et al., 2001; Merritt et al., 2001), NGC205 (Valluri et al., 2004), and M32 (Verolme et al., 2002). Some galaxies exhibit a compact stellar nucleus (with half-light radius r _h ∼ 2 − 4 pc) in their center. This is ∼ 20 times brighter than a typical globular cluster (Kormendy and McClure, 1993; Butler and Martínez-Delgado, 2005). In the Virgo and Fornax Clusters ∼ 25% of dE galaxies contain such nuclei (Binggeli et al., 1985; Ferguson, 1989; Binggeli and Cameron, 1991), but the observations with the Hubble Space Telescope revealed that these structures are far more common: about 50–80% of the less luminous galaxies contain a distinct nuclear star cluster (Carollo et al., 1998; Matthews et al., 1999; Böker et al., 2002; Balcells et al., 2003; Graham and Guzmn, 2003; Lotz et al., 2004; Grant et al., 2005; Côté et al., 2006). Nuclear star clusters are not a replacement for black holes. On the contrary low mass galaxies (10⁹ − 10¹⁰ M_⊙) show a high incidence of nuclear star clusters coexisting with massive black holes (Greene et al., 2020). However, nuclear star clusters are rare in high mass galaxies (Graham and Spitler, 2009), suggesting that the growth of BH during activity may lead to the demise of the star cluster itself (Antonini et al., 2019).

The low mass end of the relation M _BH—σ _⋆ for quiescent galaxies is still poorly sampled and of uncertain interpretation (Graham and Scott, 2015). Tidal disruption events (TDEs) provide an independent method for M _BH estimation. First, TDEs are luminous flares predominantly detected in quiescent galaxies (very few events have been detected in AGN, as the luminosity of the AGN obliterates the brightness increase associated with the TDE). Flares are produced by the tidal debris that fall back toward the black hole and that form an accretion ring or a disk around an otherwise inactive black hole. Second, a TDE can take only with relatively low black hole masses, M _BH ≲ 10⁸ M_⊙ for a solar-mass star (Gezari, 2021) to avoid that the star crosses the black hole event horizon without being tidally disrupted. The central BH mass is recovered via synthetic multi-band optical light curves based on hydrodynamical simulations of polytropic tidally disrupted stars (Mockler et al., 2019). The method is not yielding yet a good agreement with other M _BH estimates, as it is not clear which parameter should be correlated with M _BH, in a model in which the TDE luminosity is powered by fall-back accretion (Gezari, 2021).

Stellar dynamical determinations for AGN have been possible only for weakly or mildly active Seyfert 1 galaxies. Presently, only ≲ 100 dynamical M _BH measurements have been obtained by modeling stellar kinematics of quiescent and active nuclei [e.g. (Saglia et al., 2016; Kormendy and Ho, 2013)]. A case in point is the intermediate Seyfert 1 galaxy NGC3327 (Davies et al., 2006) which illustrates the complexity of the nuclear regions of a mildly active AGN, even if an application of the Schwarzschild method of orbit superposition allowed for a meaningful estimate of the stellar dynamic mass ∼ 10⁷ M_⊙. A Population B source (see §14.2 and Figure 12) radiating at modest Eddington ratio, NGC3227, shows a nuclear stellar distribution within a few parsecs of the central black hole affected by intense bursts of star formation occurring in its recent past. Similar considerations apply to the stellar dynamics results on MCG–6–30–15 (Raimundo et al., 2013). In general, stellar populations in the closeness of the active nucleus are not easy to model, also because of the uncertain distribution of obscuring dust, even in the least active nuclei, such as NGC4258. Nonetheless the stellar dynamical and maser masses agree very well for this source (Siopis et al., 2009). A dynamical mass estimate for Cen A also agrees with the estimate derived from cold molecular H₂ gas (Neumayer et al., 2007; Cappellari et al., 2009), suggesting that molecular gas could provide mass estimations as accurate as the ones based on stellar dynamics.

The most reliable method to probe distances within r _h from the black hole of galactic nuclei that are currently mildly active involves observation of H₂O masers (Miyoshi et al., 1995; Herrnstein et al., 2005; Greene et al., 2010). The H₂O emission profile shows a radial velocity “cusp” at distances where the velocity field is governed by the gravity of the black hole i.e., when r ≲ r _h. This method is not exempt by issues that could bias the results. A maser disk with Keplerian rotation could have a non-negligible disk mass comparable to the black hole mass (Huré et al., 2011), and the effects of disk self-gravity might lead to large systematic errors in the derivation of the black hole mass (Kuo et al., 2018). Disturbed morphology and kinematics are an ubiquitous feature of maser systems especially above Eddington ratio λ _Edd ≳ 0.1 (Kuo et al., 2020).

Outside of the local Universe (i.e., beyond 2.5 Mpc), VLBI observations of mega-maser systems can probe within the sphere of influence for BHs down to 10⁶ M_⊙ even at 100 Mpc (van den Bosch et al., 2016). Mega-masers are more frequently associated with active galaxies (Constantin, 2012), and they should be more common at high redshift (van den Bosch et al., 2016). The exploitation of mega-masers is however difficult, because the maser signal of high-redshift sources is lost in noise, and major surveys have until now failed to detect a mega-maser in the wide majority of galaxies. Mass determinations based on the resolved maser systems are completely independent of any other method, and best suited to cross-check the estimates obtained from stellar and gas dynamics. Several H₂O masers have been detected from the circumnuclear regions of quasars also at relatively high redshift (Broome et al., 2019; Stacey et al., 2020). The hope to go beyond modest distances rests with SKA—because of its unprecedented μJy sensitivity—and in future space based radio interferometry with μarcsec spatial resolution. Assuming Earth-space baselines of about 30,000 km, angular resolution of 2 μ-arcsec would be achievable at 8 GHz (Taylor et al., 2007). This angular scale corresponds to a projected linear size of ≈ 0.02 pc at z = 2, therefore allowing to probe within the BH sphere of influence even at the remote epochs of the “cosmic noon,” a key epoch of galaxies AGN when a population of most luminous quasars was shining bright and producing the maximum feedback on their host galaxies.

ALMA is today the most powerful tool to yield M _BH also for AGN [e.g. (Onishi et al., 2015)]. However, the CO J = 2-1 kinematics in a sample of nearby AGN reveals noncircular motions in the inner kiloparsec of all galaxies in the sample, although molecular gas and stellar kinematics show an overall agreement. The CO observations of nearby radio galaxies detect molecular disks, but also caution about the possibility of asymmetries and disruptions due to interactions with the radio jet (Ruffa et al., 2019).

Several studies have employed the capabilities of STIS on board HST to study the dynamics of line-emitting gas in proximity of the central black hole (Pastorini et al., 2007), or the sub-arcsec spatial resolution of imaging spectrometers or IFU units operating with adaptive optics (Hicks and Malkan, 2008). The concern is that radiation forces within the inner 100–1,000 pc of the central black hole may be affecting the dynamics of the line-emitting gas even more than in the case of cold gas dynamics, especially if the AGN is radiating at high Eddington ratio (Zamanov et al., 2002; Xu and Komossa, 2010; Marziani et al., 2016a; Berton et al., 2016; Schmidt et al., 2018).

Spectro-astrometry is another promising tool: the approach is based on the different photocenter positions of emission lines at different velocities (Gnerucci et al., 2011). Although a relatively modest spectral resolution ( ∼ 10 km s⁻¹) is sufficient, sub-arcsec spatial resolution is required, obviously the higher the better, achievable only from space or from ground using active optics (Abuter et al., 2021). This is the approach exploited by GRAVITY, an instrument of the Very Large Telescope Interferometer (VLTI) (Abuter et al., 2017). After first light in 2017, GRAVITY detected a spatial offset (with a resolution of 10 micro-arcseconds corresponding to approximately 0.03 parsecs) between the red and blue centers of the Paschen-α line of 3C273 (Sturm et al., 2018). This offset corresponds to a gradient in velocity and implies that the gas is orbiting the central supermassive BH. With the new capabilities of GRAVITY (Abuter et al., 2017) and with the wave-front corrections of an adaptive optics system, it will be possible to repeat this feat in many low-z type-1 (i.e., broad line) AGN (Bosco et al., 2021). The broad line region velocity field has been spatially resolved and modeled even in NGC3783 (Amorim et al., 2021) to provide an M _BH estimate, although this achievement will likely remain restricted to low-z Seyfert-1 nuclei for the time being.

If we exclude masers, for which BH masses can be inferred from rotation (Ferrarese and Ford, 2005) and spectro-astrometry, for the vast majority of Type-1 AGN the BH masses are derived from the (presumed) virial motions of the broad line region (BLR) gas clouds orbiting in the vicinity of the central compact object. If the motion in the emitting gas is in virial equilibrium, we can write the central black hole mass M _BH as: M BH = r BLR δ v K 2 G . (35) Here δv _K is the virial velocity module, r _BLR the radius of the BLR, G the gravitational constant. Eq. 35 can be useful if we can relate δv _K to the observed velocity dispersion, represented here either by the dispersion σ or by the FWHM of a suitable broad emission line: M BH = f S r BLR F W H M 2 G (36) via the structure factor (a.k.a. form or virial factor) whose definition is given by: δ v K 2 = f S F W H M 2 . (37)

Mildly ionized gas dynamics i.e., gas motions within the broad line regions of type-1 AGN, is the basis of the estimate of the M _BH for large samples of quasars up to the highest z, following Eq. 36. In addition to a measure of the virial broadening provided by the emission line width, a measure of the line-emitting gas distance from the central black hole is needed. Under the assumption that the main source of line emission is provided by photoionization (Shuder, 1981), the distance is measured by the time lag of the emission lines with respect to continuum variations (Peterson, 1993): r _BLR ≈ cτ, where τ is the time delay. Recent observations measure τ as a function of wavelength across the line profile in an attempt to resolve the velocity field of the emitting region (Brotherton et al., 2020; Williams et al., 2021). The “reverberation mapping” technique has been described in several reviews that also include a critical discussion of the technique shortcomings (Horne et al., 2004; Marziani et al., 2006; Peterson, 2014). The r _BLR estimates have been carried out mainly for the HI Balmer line Hβ for ∼ 100 type-1 AGN, recently supplemented by the monitoring of the SDSS Reverberation Mapping Project that yielded data for 144 quasars (Li et al., 2017). The reverberation mapping determinations of r _BLR offer a sort of primary step over which a correlation between r _BLR and luminosity is built (Section 13.4), in turn instrumental to the determination of the M _BH in large samples of quasars (Section 13.5).

A correlation between the radius of the emitting regions and continuum luminosity is expected on the basis of the spectral similarity of quasars. Even if this is an oversimplification, we observe always the same lines, and their relative intensities change only within a limited range, also in response to continuum variation. The ionization parameter should remain roughly constant, implying that r ∝ L ^a , with an exponent a at any rate close to 0.5 (Kaspi et al., 2000; Bentz et al., 2013). The scaling relation has been derived from spectroscopic monitoring of emission lines (mostly the HI Balmer line H β) that yield the time delay τ of the emission line response to continuum variations (Peterson, 1993; Peterson, 2017). A sufficient number of sources is available for a correlation analysis since the early 2000s (Kaspi et al., 2000). The consideration of various aspects (host galaxy subtraction and removal of the line narrow component believed to be emitted in a different region) and the increase of the number of monitored sources has led to a standard r − L relation with an exponent consistent with 0.5 within the uncertainties (Bentz et al., 2006; Bentz et al., 2009; Bentz et al., 2010; Bentz et al., 2013). However, the r − L relation suffers of significant scatter because it was derived neglecting the diversity of type-1 quasars organized by the quasar Main Sequence (MS, see below), and is biased in favor of sources radiating at a relatively low Eddington ratio. It is not difficult to account for this preferential selection: such sources are relatively low accretors and therefore more prone to variability associated with an unsteady accretion flow. Recent work (Du et al., 2016a; Du and Wang, 2019) has shown that sources that radiate at high L/L _Edd significantly deviate from the correlation of (Bentz et al., 2013): their radius is shorter than the one expected on the basis of their luminosity. Including high Eddington ratio sources in the correlation creates a cluster of data points that increases the scatter in the correlation. Figure 10 shows that the relation for sources radiating at the highest value of the Eddington ratio is significantly offset from the one of other spectral types along the quasar main sequence discussed in Section 14.2. Linear combination with the dimensionless accretion rate (i.e., the mass accretion rate normalized by the Eddington accretion rate) or Eddington ratio leads to a significant reduction of the scatter (Du and Wang, 2019; Martínez-Aldama et al., 2020).

The virial theorem can be conveniently rewritten as log  M _BH = α  log  L+β log FWHM + γ, where β = 2. The luminosity term comes from the use of the radius—luminosity relation, r − L ^a . Several different scaling laws based on this expression have been defined for the width of different lines, and for different continuum and line luminosity as well. The most widely used has been perhaps the one formulated by Vestergaard and Peterson (2006) for H β and continuum luminosity at 5,100 Å.

The main underlying assumptions in the use of the virial theorem are that the broadening is due to Doppler effect because of the line-emitting gas, and that the velocity field is such that the emitting gas remains gravitationally bound to the black hole. Early UV and optical inter-line shift analysis provided evidence that not all the line-emitting gas is bound to the black hole (Gaskell, 1982; Tytler and Fan, 1992; Brotherton et al., 1994; Marziani et al., 1996; Leighly and Moore, 2004). The emerging scenario is that outflows are ubiquitous in AGN, they occur under a wide range of physical conditions, and are detected in almost every band of the electromagnetic spectrum and on a wide range of spatial scales, from few gravitational radii to tens of kpc (e.g., (Capetti et al., 1996; Colbert et al., 1998; Everett, 2007; Carniani et al., 2015; Cresci et al., 2015; Bischetti et al., 2017; Komossa et al., 2018).

For z > 4, M _BH estimates historically rely on the Civ λ1549 high-ionization line, and the highest-z sources appear almost always high-accretors (Bañados et al., 2018; Nardini et al., 2019). The source of concern is that high-ionization lines such as Civ λ1549 are subject to a considerable broadening and blueshifts associated with outflow motions already at low redshift (Coatman et al., 2016; Sulentic et al., 2017; Marinello et al., 2020a; Marinello et al., 2020b). Overestimates of the virial broadening by a factor as large as 5–10 (Netzer et al., 2007; Sulentic et al., 2007; Mejía-Restrepo et al., 2016; Mejía-Restrepo et al., 2018) for SMBHs at high z may pose a spurious challenge to concordance cosmology (Trakhtenbrot et al., 2015) and lead to erroneous inferences on the properties of the seed BHs believed to be fledgling precursors of massive BHs. The solution is either to carry out H β observations at high redshift (a feat that is becoming easier as more NIR spectrometers are being installed at the focus of large telescopes) or to use a surrogate line whose profile is also virially broadened. The Al iii λ1860 and Ciii] λ1909 lines could be much robust estimator of the BH mass. These lines, in a blend at 1900 Å, can be easily observed with optical spectrometers up to redshift z ≈ 4. Similar considerations apply to the use of Mgii λ2800 (Shen and Liu, 2012; Trakhtenbrot and Netzer, 2012) which however can be observed only up to z ≈ 2.5 without the use of NIR spectrometers. Another approach has been to apply corrections to the Civ λ1549 line width (Coatman et al., 2017; Marziani et al., 2019), although such corrections, to be effective, require the knowledge of the quasar rest frame, which remains poorly known from rest-frame UV observations only. Shen and Liu (2012) propose scaling laws in which the virial assumption is released that is, with β ≠ 2. For Civ λ1549, this means correcting for effect associated with the emission component due to an outflow, which overbroadens (and shifts) the line. The scaling law introduced by Park et al. (2013) follows this approach assuming γ = 0.5, that is, a FWHM dependence that is much weaker than the one of the virial law. The scaling law suggested by Park et al. (2013) applied to a high luminosity sample properly corrects for the overbroadening of the Civ λ1549 line profiles of high Eddington ratio sources of a high-luminosity sample, but overcorrects the width in case of sources radiating at modest Eddington ratios, yielding a large deviation from the H β-derived M _BH values (Marziani et al., 2019).

There is a general consensus that most galaxies host massive BHs that went through phases of activity. This latter had a role in the BH growth and in the regulation of the SF activity of the host galaxy by means of wind/jet driven feedback mechanisms (Shankar et al., 2009a; Shankar et al., 2009b; Alexander and Hickox, 2012). The theoretical models show that an AGN and its host may coevolve (Silk and Rees, 1998; Granato et al., 2004), leading to characteristics (such as the M _s,sph /M _s ratio and/or the central stellar velocity dispersion σ) related to black hole mass (M _BH ).

An early answer to the question “do galaxies hosting an AGN share the same M _BH − M _Bulge correlation of normal galaxies?” was affirmative: AGN have the same BH-bulge relation as ordinary (inactive) galaxies (Wandel, 2002). Fast forward 20 years, there is not yet an established view. A most recent work, based on state-of-the-art surface photometry, and spatially resolved kinematics to measure σ _⋆, find that correlations between M _BH and host galaxy properties hold for AGN within the limits of an intrinsic scatter 0.2–0.4 dex, and are consistent with the ones of quiescent galaxies (Bennert et al., 2021).

Recent works also point toward a complex scenario involving selection biases (Schulze and Wisotzki, 2011) and a better appreciation of the active galaxies diversity. We may represent the distribution of objects in the M _BH—M _bulge (or σ _⋆) diagram by the bivariate distribution function of bulge mass M _bulge and M _BH Ψ(M _BH,M _bulge). The Ψ distribution can be factorized as Ψ = γ(M _BH | M _bulge)ϕ(M _bulge) where ϕ(M _bulge) is the spheroid mass function and γ represents the M _BH—M _bulge correlation i.e., the probability of having the black hole mass M _BH for a given M _bulge. A correct evaluation of γ(M _BH | M _bulge) relies on: 1) the knowledge of ϕ(M _bulge), which is not a trivial task to achieve even in the local Universe, and needs a separate consideration of purely spheroidal (i.e., diskless) galaxies and galaxies with pseudo-bulges or with a bulge/disk system; 2) the absence of biases affecting γ(M _BH | M _bulge).

Both the determinations of M _BH and bulge parameters are challenging, when derived from conventional optical and NIR measurements. At present, black hole masses for type-1 AGN are more frequently derived through the so-called, single epoch virial broadening estimation i.e., through the measurement of the radial velocity broadening term that appears squared in Eq. 36 from single epoch spectra. In practice, it is the measurement of the FWHM or σ ⁶ of broad emission lines [e.g. (McLure and Jarvis, 2002; Vestergaard and Peterson, 2006)]. To obtain M _BH, an estimate of the radius r _BLR is also needed, and a rather poorly defined scaling law of r _BLR vs luminosity is applied.

The bulge estimates in AGN samples are hampered by the luminous source associated with the active nucleus, which may well outshine the entire galaxy. The M _bulge has been frequently computed from the host galaxy luminosity (Peng et al., 2006a; Peng et al., 2006b; Treu et al., 2007; Schramm et al., 2008; McLeod and Bechtold, 2009; Bennert et al., 2010; Decarli et al., 2010; Targett et al., 2011). However, type-1 AGN remain offset from inactive galaxies in the M _BH—L _bulge relation: AGN have more luminous bulges at a given black hole mass (Nelson et al., 2004; Bennert et al., 2021). There are evidences that the M _BH − L _bulge relation defined by quiescent BH samples differs from that defined by the galaxies in the SDSS (Bernardi et al., 2007). Interestingly, the offset is larger for AGN of larger Eddington ratio (Barth et al., 2021). This suggests that the central regions of galaxies hosting an AGN have, in general, lower mass-to-light ratios than inactive galaxies, most likely for the presence of a young stellar population in the bulge of active systems (Kim and Ho, 2019).

The σ _⋆ has been measured either directly i.e., from the width of absorption lines associated with the stellar component of the host galaxies (Woo et al., 2006; Woo et al., 2008b; Shen et al., 2008) or by using the widths of narrow emission lines a proxy of the stellar velocity dispersion (Shields et al., 2003; Shields et al., 2006; Salviander et al., 2007). The latter approach is fraught from systematic effects. In the case of quasars and AGN radiating at moderate and high Eddington ratios, the [Oiii] λ5007 broadening is strongly affected by non-virial motions (Zamanov et al., 2002; Marziani et al., 2006; Mathur et al., 2011; Marziani and Sulentic, 2012; Cracco et al., 2016).

More reliable results for dynamical mass measurements of the host galaxy from spatially resolved images have been obtained with adaptive optics (Inskip et al., 2011). CO emission profiles have been used to estimate dynamical masses for individual objects since the early 2000s (Walter et al., 2003) even at fairly high redshift, and nowadays ALMA is rapidly adding to the available dynamical mass measurements for the host galaxy [e.g. (Tan et al., 2019; Molina et al., 2021)], considering that the velocity field of the molecular gas is often regular and consistent with rotation. State-of-the-art surface photometry of the AGN host galaxies in the NIR achieves decomposition in spheroid, disk, and bar component, as most of the host of nearby Seyfert galaxies are of morphological type Sa/SBa. As mentioned, a most recent work did not detect significant differences in the scaling with M _BH and σ _⋆ between active and nonactive galactic nuclei (Caglar et al., 2020; Bennert et al., 2021), and did not find difference between pseudo and classical bulges or barred and nonbarred galaxies in the M _BH—M _bulge relation (Bennert et al., 2021), although this result is still controversial (Kormendy et al., 2011; Ho and Kim, 2014). In addition, γ is still computed with the single epoch technique, without consideration of the diversity in accretion structure (and hence virial factor) that is expected in type-1 AGN samples. For type-1 active nuclei radiating at Eddington ratio above 0.01, the geometry and structure of the emitting region are affected by the accretion mode, which in turn affects the expression of the virial factor that is dependent on kinematics, geometry, and viewing angle (Collin et al., 2006; Park et al., 2012; Mejía-Restrepo et al., 2016; Mejía-Restrepo et al., 2018; Shankar et al., 2019). A study separately considering sources in different accretion modes and the statistical bias introduced by orientation effects is not available as yet.

Keeping the attention focused on γ, the radius of influence r _h is of the order of parsecs, and insufficient resolution may prevent reliable BH mass estimates or forces to target only the largest BHs (Gültekin et al., 2009; Gültekin et al., 2011), leading to a selection effect that yields an increase in the M _BH—σ relation for quiescent galaxies by a factor of a few (Shankar et al., 2016; Shankar et al., 2019). AGN will on average host more massive BHs than in the volume-limited case (Lauer et al., 2007), determining a Malmquist bias toward more massive BHs at a given spheroid mass, shifting γ upward and causing an offset in the zero-point of the M _BH—M _bulge relation. There are competing effects: the fraction of active galaxies among SMBHs varies considerably with mass (high-mass BHs are likely less active than low-mass BHs (Schulze and Wisotzki, 2011)). The strength of the bias depends on the limit in luminosity, the shape of the distribution function of spheroids, the scatter of the M _BH-M _bulge relation, and the Eddington ratios. If, as mentioned, the active fraction decreases as the BH mass increases, then for a given spheroid mass it will be more probable to find small-masses BHs in an AGN sample, causing a bias toward lower M _BH/M _bulge ratios, and a change in the slope of the relation.

The low-M _BH end of the correlation is especially problematic, as it is for quiescent galaxies. Narrow-line Seyfert 1s nuclei (NLSy1s, low-z type-1 AGN several of which are accreting at high rate (Marziani and Sulentic, 2014a)), often hosted in dwarf high surface brightness galaxies (Krongold et al., 2001) and in barred spirals (Crenshaw et al., 2003; Ohta et al., 2007), possess under-massive BHs (Mathur et al., 2001; Chao et al., 2008). NLSy1 nuclei often reside in disk-dominated galaxies with pseudo-bulges (Orban de Xivry et al., 2011; Mathur et al., 2012; Ermash and Komberg, 2013; Olguín-Iglesias et al., 2017; Järvelä et al., 2018; Doi et al., 2020). These types of bulges are more closely associated with the evolution of disks and may be typical of systems that did not experience a minor or major merger capable of leading to a real bulge development. Several studies found that disk-dominated galaxies deviate from the M _BH—M _bulge correlation, and, if considered as a distinct class, may not follow a M _BH—M _bulge correlation (Kormendy et al., 2011; Davis et al., 2018; Sahu et al., 2020). However, if one applies a correction for the disk component, and considers only the bulge, the AGN in the low black hole mass ranges M _BH ≲ 10⁸ M_⊙ might follow a relation consistent with the local M _BH—M _bulge correlation (Bennert et al., 2011; Sanghvi et al., 2014). At any rate, the relation between M _BH and σ _⋆ or M _bulge should be taken with special care in particular in the lower M _BH range. Relatively few objects are obscured type-1 AGN. Chandra observations are detecting a wealth of black holes in star-forming galaxies, in the range between 10⁶—10⁷ M_⊙, even at high z (Mezcua et al., 2016; Fornasini et al., 2018; Zou et al., 2020). They are low mass by supermassive black hole mass standards, and most likely still growing in an obscure phase. It is not known how they would be located in the M _BH—M _bulge plane. These elusive AGN are potential targets for JWST (Satyapal et al., 2021).

At the time of its discovery, the luminous quasar HE0450-2,958 appeared as an oddity: a quasars without a host galaxy! (Magain et al., 2005). Understandably enough, the source attracted a lot of interest, and perhaps even a revival of the noncosmological interpretation of quasar redshifts (Arp et al., 1979; Sulentic and Arp, 1983; Sulentic and Arp, 1987). HE0450-2,958 appears hosted by a galaxy much fainter than that inferred from the correlation between BH mass and bulge luminosity (Kim et al., 2007). In the case of quiescent galaxies, compact dwarf galaxies whose BH has a mass reaching even 15% of the total galaxy mass (Reines et al., 2014; Seth et al., 2014; van Loon and Sansom, 2015) are observed. A possible explanation is that their outer parts may have been stripped by repeated encounters with other galaxies and produced an ultra-compact dwarf galaxy. The EAGLE cosmological and hydrodynamical simulations suggest that these kinds of objects are outliers resulting from the combination of stellar tidal stripping and the early formation epoch, which leaded to a rapid BH growth at high redshift, with the first mechanism being the most relevant for 2/3 of these sources (Barber et al., 2016). However, the disk/bulge decomposition is a delicate procedure. A careful reanalysis of the most striking cases, Mrk1216, NGC1277, NGC1271, and NGC1332, suggests that a proper reevaluation of the disk size with an ensuing increase in spheroid mass will bring these sources in better agreement with the M _BH—M _bulge relation (Savorgnan and Graham, 2016). The case of HE0450-2,958 has not been fully explained to date. Past works have considered intriguing lines of evidence suggesting high L/L _Edd and BAL outflow (Merritt et al., 2006; Lipari et al., 2007). However HE0450-2,958, which appears as a mini-BAL from a FOS spectrum, shows modest optical Feii emission, and a spectrum similar to the one of PG1211 + 143 (Merritt et al., 2006). According to the main sequence trends (see §14.2), the object should not be highly accreting (Marziani and Sulentic, 2014a; Du et al., 2016b). It is also unlikely that HE0450-2,958 is a recoiling black hole ejected by a companion galaxy at approximately 7 kpc of projected linear distance, on the ground of the strong narrow line emission of [Oiii] λλ4959,5007 (Merritt et al., 2006). HE0450-2,958 does not appear as an extraordinary powerful quasar. The upper limits on the host galaxy luminosity are not very constraining, so that this object could be well within the limits set by the scatter in the M _bulgevM _BH correlation (Kim et al., 2007).

However, recoiling black holes—provided that they are the active member of the binary, as suggested by numerical simulations (Nguyen and Bogdanović, 2016)—may systematically lower M _BH and ultimately increase the scatter of the observed BH–host galaxy bulge relation due to ejected BHs (Volonteri, 2007). Recoiling BHs have lower masses than their stationary counterparts, but the deficit in mass depends on kick speed and merger remnant properties (Blecha et al., 2011). The effect is of an overall downward shift in the normalization and an increase of the scatter in the M _BH—M _bulge relation: the offset between the stationary and the recoiling BH population can reach δ  log  g ≈ 0.4 dex, if the rotational velocity of the secondary BH is close to its escape velocity. The amplitude of the downward offset depends on the recoil velocity as well as on the accretion history of the stationary black hole, and can be lower, yielding a δ  log  ϕ ≈ 0.2 dex. This scenario is not as yet contextualized: a large fraction of type-1 AGN shows evidence that they do not host a sub-parsec binary black hole with a significant mass ratio between the secondary and the primary (say q ≳ 0.1). Conclusive evidences in favor of such binary systems are very rare at the time of writing.

Active galactic nuclei and quiescent bulge-dominant galaxies do not show strong evidence of evolution in the M _BH—M _bulge relation up to z ∼ 0.6 − 1 (Salviander et al., 2007; Schulze and Wisotzki, 2014; Li et al., 2021). At higher redshift, there is an increasing evidence of evolution, in the sense of high-z SMBHs that are overmassive at a given bulge mass than expected from the local scaling relation (McLure et al., 2006; Decarli et al., 2018). Between redshifts 1 and 2, Merloni et al. (2010) suggested a significant increase of the M _BH/M _Bulge ratio ( ∝ (1 + z)^0.68). Studies at even higher redshift used the velocity dispersion of the gas as a proxy of the stellar velocity dispersion and dynamical mass measurement from inclined disk models (Vayner et al., 2021). They suggest over-massive black holes (Targett et al., 2011) with respect to the local scaling law. The most recent results confirm that quasars host galaxies are under massive relative to M _BH, and detect a large difference, even by an order of magnitude, with systems at redshift in between 1.4 and 2.6 residing off the local scaling relation. Several quasar host galaxies have been resolved in their [C II] emission on a few kpc scale at redshift ≈6. Even in this case, the dynamical mass estimates for the host galaxies give masses more than an order of magnitude below the values expected from the local scaling relation (Decarli et al., 2018), in agreement with the results for galaxies at z ≈ 7 derived from cosmological hydrodynamical simulations (Marshall et al., 2020).

The evolution of the M _BH—M _bulge relation with the cosmic epochs can be interpreted in several ways: the most straightforward is that of a rapid growth of SMBHs at high redshift (Lupi et al., 2021). Also a variation of structural properties of AGN hosts remains possible (Shankar et al., 2013b; Zhu et al., 2021): elliptical galaxies are not really monolithic spheroids, but have undergone significant late-time dissipation-less assembly. There are intriguing caveats with the interpretation of a rapid black hole growth. First, very massive seed black holes need to be formed at z ≈ 20 to account for masses ∼ 10⁹ M_⊙ observed at redshift z ≳ 4 ((Volonteri, 2010; Trakhtenbrot, 2020), and references therein). Second, BH masses (unlike the masses of galaxies!) can only increase with cosmic epoch. If the merger-driven hierarchical scenario that implies the parallel growth of bulges and BHs is taken literally, the larger M _BH/M _Bulge ratio at high z means that mergers affect more bulge than BH masses (at cosmic epochs associated with z ≳ 1), an implication consistent with the anti-hierarchical growth and downsizing of the nuclear activity at low-z (Hirschmann et al., 2012). If pseudo-bulges follow the same M _BH scaling relations as that of classical bulges [e.g., Bennert et al., 2021], hierarchical growth might not be the only mechanism that drives the relation between M _BH—M _bulge: in spiral galaxies, secular evolution might lead to a parallel growth of bulge and central black hole. Clearly, this issue should be analyzed in connection to ongoing star formation properties of the pseudo-bulge hosts (Zhao et al., 2021). Host and black hole properties are different for different masses, and the relation between galaxy color and black hole mass is different for the red and blue sequence quiescent galaxies, suggesting different channels of black hole growth for the two sequences (Dullo et al., 2020).

In conclusion, AGN with “coreless” elliptical/bulge-dominated hosts may straightforwardly follow a relation similar to the one of normal galaxies. In other words, the M _BH—M _bulge relation may strictly hold for massive evolved systems, also if the nucleus is active, in a form that is as yet indistinguishable from the one of quiescent galaxies. It remains to be tested whether these sources could be mainly AGN accreting at relatively low rate and radiating at modest Eddington ratios (Population B, Section 14.2). Significant deviations may be associated with disk dominance, but a careful assessment of the relative disk and bulge contribution might bring the system with the over-massive BHs in agreement with the established relation (Caglar et al., 2020; Zhao et al., 2021). The local NLSy1s—all of which are Population A (Section 14.2), with a significant fraction of high accretors—are instead believed to be with black holes under massive with respect to their host masses. In this respect they are different from the high-z quasars with over-massive black holes. However, the observational properties of low-z AGN accreting at relatively high rate can still be regarded as typical of very high z quasars, when massive bulges were not yet formed, as originally suggested by Mathur (2000), Sulentic et al. (2000a). The analogy is based on the optical, UV, and X-ray AGN spectroscopic properties that are mainly governed by the Eddington ratio. In addition, modest masses of low-z quasars can grow by a factor ∼ 10 on time scales shorter than timescale of the cosmic evolution of quasar accretion rates, and therefore bring under massive BHs in line with the M _BH—M _bulge relation (Fraix-Burnet et al., 2017).

The fundamental plane of black hole activity can be written as a correlation between black hole mass, X-ray, and radio luminosity. The correlation defines a plane in the space of parameters defined by the mass and the radio and X-ray luminosities. In its original formulation, the fundamental plane was written as (Merloni et al., 2003): l o g L R = ( 0.60 ± 0.11 ) log L X + ( 0.7 8 − 0.09 + 0.11 ) log M B H + 7.3 3 − 4.07 + 4.05 . (38)

The scatter is large, implying that a fourth variable might be involved, for instance black hole spin (Ünal and Loeb, 2020). The salient point is however that the relation holds over a huge range of black hole masses, from a few times solar (i.e., from the domain of the so-called micro-quasars) to the largest black hole masses detected in the Universe ∼ 10¹⁰ M_⊚ [e.g. (Schindler et al., 2021; Valtonen et al., 2012)]. ⁷ It is remarkable that also stellar-mass black holes exhibit relativistic jets, as spectacularly demonstrated by the relatively recent discovery of superluminal motion in a Galactic black hole candidate by Mirabel and Rodríguez (1994). The self-similarity expressed in Eq. 38 allows for an invariant jet model and a simple relation between M _BH and radio power (Heinz and Sunyaev, 2003). The self-similarity notwithstanding, there is a nonlinear relation between BH mass and radio power, with P _ν ∝ M BH 1.3 − 1.4 , implying that the radio emission normalized to the bolometric luminosity should be much higher for AGN than for microquasars. In the framework of the model of Heinz and Sunyaev (2003), flat spectrum synchrotron jet emission is produced by an inefficient accretion mode. The fundamental plane of black hole activity refers to sources accreting at very low rate (dimensionless accretion rate m ̇ ≲ 0.01), and radiating below a few hundredths of their Eddington luminosity (Gültekin et al., 2019). This means that the relation is best suited for sources such as micro-quasars (i.e., stellar-mass black hole candidates in the low state (Mirabel and Rodríguez, 1994)) and BL Lac objects, in which both radio and X-ray emissions are ultimately associated with the relativistic jet. ⁸

The quasar main sequence is, in many ways, analogous to the FP for black holes in a different accretion mode sustained by higher m ̇ ∼ 0.01 − 1 . The formulation is rather different, and follows a different discovery path based on the statistical analysis of sources that are predominantly radio-quiet. The quasar main sequence (MS) is defined from the first Eigenvector (E1) that was originally identified by a PCA of about 80 Palomar-Green (PG) quasars and associated with an anti-correlation between the strength of optical Feii emission measured from the prominence of the emission blend centered at λ 4,570 Å (Fe ii λ4570) with respect to H β (R _FeII = I (Fe ii λ4570)/I(H β)) and FWHM of H β (Boroson and Green, 1992). The E1 MS has withstood the test of time [Figure 11 (Sulentic et al., 2000b; Zamfir et al., 2010; Popović and Kovačević, 2011; Kruczek et al., 2011; Grupe and Nousek, 2015; Shen and Ho, 2014),], and the main optical trend shown in Figure 11 has been confirmed by samples of more than two order of magnitude larger in size than the original one (Shen and Ho, 2014). The importance of Feii stems from its extensive emission from UV to the IR that can dominate the thermal balance of the low-ionization BLR. The FWHM(H β) is associated with the velocity field in the low-ionization BLR, most likely predominantly virialized (Peterson and Wandel, 1999). These two parameters are related to the physical conditions and to the dynamics of the emitting regions, which are in turn influenced by the accretion mode of the central black hole, and its evolutionary stage.

Trends associated with the MS have been extended to the radio (Sulentic et al., 2003; Shen and Ho, 2014; Zamfir et al., 2008; Ganci et al., 2019), FIR (Wang et al., 2006; Ganci et al., 2019), IR (Dultzin-Hacyan et al., 1999; Loli Martínez-Aldama et al., 2015; Panda et al., 2020), UV (Sulentic et al., 2000c; Reichard et al., 2003; Bachev et al., 2004; Sulentic et al., 2006; Śniegowska et al., 2020; Baskin and Laor, 2005; Richards et al., 2002; Richards et al., 2005) and X-ray domain (Wang et al., 1996; Grupe et al., 2001; Bensch et al., 2015), and to optical variability as well (Mao et al., 2009; Bon et al., 2018). Table 1 of Fraix-Burnet et al. (2017) provides a detailed list of the various parameters that have been measured in the various frequency domains. A summary description of the trends and a justification for the two quasar populations are also provided by several authors (Sulentic et al., 2008; Sulentic et al., 2011; Sulentic and Marziani, 2015). A nonlinear decay curve provides a quantitative description of the main sequence in the FWHM—R _FeII plane (Wildy et al., 2019).

The distribution of the data in the plane R _FeII—FWHM(H β) makes it expedient to define spectral types [Figure 12 (Sulentic et al., 2002; Shen and Ho, 2014)]. This provides the considerable advantage that a composite spectrum within each bin could be representative of objects in similar physical conditions. In alternative, a prototype object can be defined for each spectral type and used to analyze systematic changes along the quasar MS. It is also expedient to distinguish between two populations: Population A made of sources with FWHM(H β) ≤ 4,000 km s⁻¹ and Population B (broader) with FWHM(H β) > 4,000 km s⁻¹. Extreme Population A are quasars with R _FeII ≳ 1 and extreme Pop. B with undetectable Feii emission and the broadest Balmer lines (extreme FWHM H β can reach ∼ 15, 000 − 20, 000 km s⁻¹ (Eracleous and Halpern, 2003; Strateva et al., 2003; Eracleous and Halpern, 2004)). Basically, Population B includes sources termed as “disk dominated,” where radiation forces exert a modest influence on the overall dynamics of the gas (Richards et al., 2002), while Population A is made of quasars radiating at relatively high Eddington ratio L/L _Edd ≳ 0.2, for which radiation forces are able to maintain a wind that leads to easily identified systematic wavelength displacements toward the blue with respect to the quasar rest frame in the high-ionization emission lines (Gaskell, 1982; Brotherton et al., 1994; Marziani et al., 1996; Richards et al., 2011; Coatman et al., 2016; Sulentic et al., 2017). The extreme of Population A identifies the class of “strong Feii emitters” (Lipari et al., 1993; Graham et al., 1996b). Feii emission overwhelming H β line emission (R _FeII ≳1) implies extreme Eddington ratio (L/L _Edd ∼ 1 (Marziani and Sulentic, 2014a)) and possibly super-Eddington accretion rate (Wang et al., 2014a; Sun and Shen, 2015; Du et al., 2016b; Du et al., 2016b; Panda et al., 2018; Panda et al., 2019).

However, along the entire MS, the BLR gas emitting the low-ionization lines belongs to predominantly virialized systems (Peterson and Wandel, 1999). The main asymmetries in the low-ionization line profiles can be explained in the context of a dynamical system whose velocity field is predominantly Keplerian. The single peaked, symmetric, and unshifted profile typical of Population A has been traditionally explained as due to an extended disk (Dumont and Collin-Souffrin, 1990), and the same explanation apparently remains valid in the case of extreme Pop. A AGN that are characterized by extreme high-ionization blueshifts (Leighly and Moore, 2004; Sulentic et al., 2007; Richards et al., 2012; Marziani et al., 2016b; Bischetti et al., 2017; Sulentic et al., 2017). The high Civ λ1549/H β intensity ratio of the blueshifted emission (Marziani et al., 2010) makes it possible that the H β profile remains almost symmetric and can be easily symmetrized by applying a small correction (Negrete et al., 2018). In general, the distinguishing feature of Pop. B sources, a redward asymmetric profile, can be explained by the sum of a disk contribution and emission from a larger distance (Bon et al., 2007; Bon et al., 2009). Reverberation mapping studies of lines from different ionic species have provided evidence of “ionization stratification” and velocity-resolved reverberation mapping of sources with asymmetric H β basically confirms the scenario of a Keplerian velocity field (Du et al., 2018b; Brotherton et al., 2020). The red-ward asymmetry has been interpreted as due to gravitational and transverse redshift (Bon et al., 2015; Punsly et al., 2020) or by gas clouds infalling toward the central black hole (Wang et al., 2017). At the extreme end of Pop. B sources, the profiles are often very broad and double peaked, accounted for by a bare Keplerian disk model with mild relativistic effects (Chen and Halpern, 1989; Strateva et al., 2003). So, all along the quasar MS the low-ionization lines (at variance with the high-ionization emission) appear to be predominantly associated with a bound, Keplerian dynamical system (Collin-Souffrin et al., 1988; Elvis, 2000).

Many studies still distinguish between the NLSy1s (FWHM H β ≲ 2000 km s⁻¹) and the rest of type-1 AGNs [e.g., Cracco et al., 2016], and consider NLSy1s an independent class. There is a general consensus that the limit at 2,000 km s⁻¹, albeit of historical importance, has no special meaning. The main reason behind extending the limit from 2,000 to 4,000 km s⁻¹ is that several properties of NLSy1s are consistent with the ones of “the rest of Population A” in the range 2,000 km s⁻¹ ≲ FWHM(H β) ≲ 4000 km s⁻¹. The change—in low redshift samples z ≲ 1—occurs around 4,000 km s⁻¹, not 2,000 km s⁻¹ (Cracco et al., 2016). On the converse Population A and B can be distinguished on the basis of the Balmer line profiles, and because of the amplitude of the systematic blueshift of the high-ionization lines with respect to the quasar rest frame. Composite Hβ profiles of spectral types along the MS are consistent with a Lorentzian for both NLSy1s and the rest of Population A. Other parameters (CIVλ1549 centroid, R _FeII) also span the same ranges in NLSy1s and the rest of Population A.

The governing accretion parameter accounting for the MS trends is most likely the Eddington ratio, which is related to the mass accretion rate by a monotonic albeit nonlinear relation (Mineshige et al., 2000; Sadowski, 2011; Sądowski et al., 2014). This explanation—originally suggested by Boroson & Green (Boroson and Green, 1992)—has also withstood the test of time (Marziani et al., 2001; Boroson, 2002; Ai et al., 2010; Zamfir et al., 2010; Xu et al., 2012; Shen and Ho, 2014; Sun and Shen, 2015; Panda et al., 2017; Panda et al., 2018), even if several key pieces needed to connect L/L _Edd to the observed parameters remain poorly understood to date. The evidence of a correlation between R _FeII and L/L _Edd is still made murky by the strong effect of orientation on the line broadening, affecting M _BH and L/L _Edd computations with both random and systematic errors (Marziani et al., 2019). Sun and Shen (2015) provided evidence of this based on the stellar velocity dispersion of the host spheroid (a proxy for M _BH) anticorrelation with R _FeII, implying that L/L _Edd increases with Feii. Recent approaches include a careful analysis of the role of metallicity and of density and ionization trends (Panda et al., 2018; Panda et al., 2019), and confirm L/L _Edd as the main physical parameter governing the MS “horizontal branch” along the R _FeII axis.

A toy scheme can explain in a qualitative way the occupation of the MS plane under the assumptions that Eddington ratio, mass, and an aspect angle θ (i.e., the angle between the line-of-sight and the accretion disk axis) are the parameters setting the location of quasar along the MS (Marziani et al., 2001; Marziani et al., 2018). If the BLR radius follows a scaling power-law with luminosity (r ∝ L ^a, Kaspi et al. (2000), Bentz et al. (2013)), under the standard virial assumption, then F W H M ∝ f S ( θ ) − 1 2 L M BH − a 2 M BH 1 − a 2 ∝ f S − 1 2 L 1 − a 2 L M BH − 1 2 . (39)

We can also write R_FeII as a function of (L/L _Edd) and θ, which needs to be established either empirically or theoretically. For illustrative purposes, we consider the “fundamental plane of accreting BHs” that relates L/L _Edd to R _FeII (Du et al., 2016b; Bon et al., 2020), ignoring other relevant factors, such as systematic differences in line shapes and in chemical composition along the MS (Panda et al., 2019; Śniegowska et al., 2020), and we assume that R _FeII depends on θ following a limb-darkening law (Marziani et al., 2001; Netzer, 2013).

As expected, the right panel of Figure 12 shows that θ predominantly affects FWHM H β and L/L _Edd predominantly (but not exclusively) affects R _FeII. Under the assumptions of the toy scheme the FWHM limit at 4,000 km s⁻¹ should include mainly sources with L/L _Edd ≳ 0.1 − 0.2. Sources at lower L/L _Edd are expected to be rare because they should be observed almost pole-on (for example, core-dominated radio-loud quasars whose viewing angle θ is relatively small (Marziani et al., 2001; Zamfir et al., 2008)), and the probability of observing a randomly oriented source at an angle θ between the symmetry axis and the line of sight is P(θ) ∝  sin θ. Even if such sources are expected to be rare, their number increases in flux-limited samples for a Malmquist bias, due to a continuum enhancement via relativistic beaming. We can say that separating Pop. A and B at 4,000 km s⁻¹ makes sense for low z samples and that, also by a fortunate occurrence, Pop. A includes mostly relatively high L/L _Edd sources.

The bolometric luminosity L can be estimated from optical or UV luminosities. ⁹ The diagram L/L _Edd vs. bolometric luminosity (Figure 13) also provides a strong rationale for the existence of two populations: only above a threshold of L/L _Edd ≈ 0.1 large shifts are observed. Data points whose high-ionization lines are strongly blue shifted with respect to the rest frame are superimposed on the distribution of Figure 13, and are clearly seen for L/L _Edd ≳ 0.1 only. This corresponds to the population A and B of Sulentic et al. (2000a), of wind and disk-dominated quasars (Richards et al., 2002), and population 1 and 2 of Collin et al. (2006). The data of Figure 13 refer to sources with large blueshift in [Oiii] λλ4959,5007, but an equivalent behavior is observed also for the blueshift of Civ λ1549. At the same time, Figure 13 (and Figure 15 as well) show the effect of a strong bias typically affecting quasar studies over a broad range of redshifts: at high z we detect only the high-luminosity sources that correspond to relatively high L/L _Edd.

Figure 11 refers to low-z (z ≲ 1) samples. A complete mapping of the MS at high L is still missing (we consider high-luminosity quasars those with bolometric log  L ≳ 47 [erg/s]): the H β spectral range is therefore accessible only with IR spectrometers to observe the H β spectral regions of high-luminosity quasars that are very rare at z ≲ 1. A significant progress is expected in the next years, since IR spectral observations covering H β of high-z and high-L quasars are becoming widespread. A systematic increase in BH mass M _BH has a corresponding increase in FWHM. If a = 0.5, the FWHM grows with M BH 0.25 , i.e. a factor of 10 for   log  L, passing from 44 (relatively low luminosity) to 48 (very luminous quasars). The trend may not be detectable in low-z flux-limited samples, but becomes appreciable if quasars over a wide interval in L are considered. At high M _BH, the MS becomes displaced toward higher FWHM values; the displacement probably accounts for the wedge-shaped appearance of the MS when large samples of quasars are considered (Shen and Ho, 2014). If we consider a limiting Eddington ratio (L/L _Edd ∼ 0.1—0.2) as a physical criterion for the distinction between Pop. A and B, then the separation based on the FWHM becomes luminosity dependent. According to the toy scheme, the FWHM of H β (or of any other virialized line) should be ∝ L / L Edd − 1 2 × L 1 − a 2 . Figure 14 shows that the ∝ L ^0.25 for the width of a low- and an intermediate ionization line. The maximum L/L _Edd should correspond to the minimum FWHM, expected to increase with luminosity as ∝ L ^0.25. If the FWHM is plotted against the luminosity, a trend-line nicely envelops the lower FWHM end of the data point distribution (Marziani et al., 2009).

Joining the fundamental plane and the main sequence trends for AGN, four main regimes can be isolated (c.f. (Giustini and Proga, 2019)) where the physics of the inner accretion and ejection is expected to change. Observationally, they range from low-luminosity AGN at extremely low accretion rates ( m ̇ ≲ 0.01 ) and Population B quasars radiating at rates 0.01 ≲ L/L _Edd ≲ 0.1 − 0.2, to Population A sources with L/L _Edd ≳ 0.1 − 0.2, and extreme Population A sources radiating close or somewhat above the Eddington limit (L/L _Edd ≳ 1). There is a close formal analogy between the FP of accreting black hole and the MS. Eq. 38 can be rewritten as an implicit relation between L/L _Edd and M _BH. Similarly the MS is a sequence in the plane FWHM H β—R _FeII that can be translated into a relation between L/L _Edd ( ∝ R _FeII) and M _BH ( ∝ FWHM). The relations of the MS are, as in the case of the FP, self-similar over 9 orders of magnitude in M _BH (Zamanov and Marziani, 2002). Radiation-driven winds appear to dominate in the high-ionization line emission in Population A and especially extreme Pop. A, reflecting the importance of the balance between radiation and gravitation forces expressed by L/L _Edd in the accretion processes of AGN (Ferland et al., 2009; Marziani et al., 2010), whereas the black hole mass is the ultimate parameter governing the energetics (Sulentic et al., 2017).

The M _BH—luminosity relation can be constructed for large quasars samples once the M _BH has been computed (Figure 15). Figure 15 shows that the distribution of quasars in the plane M _BH—L is constrained within two well-defined diagonal lines, corresponding to the L/L _Edd ≈ 0.01 and L/L _Edd ≈ 1. The empty area at the top left corner is due to inefficient radiators accreting at very low rate (Narayan and Yi, 1995), which are most often not type-1 quasars and are difficult to detect; the bottom right area is associated with sources that should be super Eddington radiators. Such sources are not expected to exist; L/L _Edd ≈ a few could be a physical limit for highly super-Eddington accretion (Mineshige et al., 2000; Sądowski et al., 2014).

Three conditions are satisfied for xA quasars: 1) constant Eddington ratio L/L _Edd, close to Eddington limit; 2) the assumption of virial motions of the low-ionization BLR, so that the black hole mass M _BH can be expressed by the virial relation (Eq. 36); 3) spectral invariance: for extreme Population A, the ionization parameter U can be written as U = Q ( H ) / 4 π r BLR 2 n H c ∝ L / r BLR 2 n H (Netzer, 2013), where Q(H) is the number of hydrogen-ionizing photons. U has to be approximately constant; otherwise, we would observe a significant change in the spectral appearance. The three constraints make it possible to derive a relation between line width (the FWHM of the H β broad component is expressed in units of 1,000 km s⁻¹) and luminosity: L ( F W H M ) = L 0 ⋅ ( F W H M ) 1000 4 e r g s − 1 (40) where L 0 depends on the square of L/L _Edd, the ionizing range of the spectral energy distribution, and a parameter directly derived from the UV spectra, the product density times ionization parameter that has been scaled to the typical value 10^9.6cm⁻³ (Padovani and Rafanelli, 1988; Matsuoka et al., 2008; Negrete et al., 2012). Until now, the FWHM of H β broad component and of Al iii λ1860 have been adopted as VBEs (Dultzin et al., 2020; Czerny et al., 2020; Marziani et al., 2020). Equation (40) implies that a simple measurement of the FWHM of a low-ionization line yields a z − independent estimate of the accretion luminosity [Marziani and Sulentic, 2014a, c.f. (Teerikorpi, 2011)].

The virial luminosity equation is conceptually equivalent to the Tully-Fisher and the early formulation of the Faber Jackson laws for ETGs (Faber and Jackson, 1976; Tully and Fisher, 1977). Recent works proposed the “virial luminosity” could provide suitable distance indicators because several emission properties appear to be extreme and stable with luminosity scaling with black hole mass at a fixed ratio (Wang et al., 2013; Wang et al., 2014a; Franca et al., 2014). The virial equation has been applied to xA quasars only (L/L _Edd ∼ 1), although it in principle could be useful for all quasars with known L/L _Edd, provided a suitable emission line broadened by virial motions is used for the luminosity computation. At present, the virial equation can be considered for all xA quasars distributed over a wide range of luminosity and redshift, where conventional cosmological distance indicators are not available (Czerny et al., 2018; Czerny et al., 2020; Marziani et al., 2020).