© The Authors 2023

1. Introduction

This study is the latest in a series aimed at demonstrating that the scaling relations (Sc–Rs) for early-type galaxies (ETGs), which are the mutual correlations between the main structural parameters of galaxies (stellar mass M_s, effective radius R_e, effective surface intensity I_e, luminosity L, and central velocity dispersion σ)¹, can be fully understood if we adopt a new perspective. The idea is that the virial theorem (VT) of the stellar systems can be coupled to the galaxy luminosity taking into account that this latter can randomly vary with time as a result of accretion–depletion events associated with mergers and/or close encounters expected in the hierarchical galaxy formation scenario in addition to the natural evolution of its stellar content. The new equation governing the luminosity is expressed by

$$\begin{array}{c}\hfill L(t)=L_0^{\prime }(t)\sigma {(t)}^{\beta }(t),\end{array}$$(1)

where L, $L_0^{\prime }$, σ, and β are all functions of time and can vary from galaxy to galaxy. This relation is formally equivalent to the Faber & Jackson relation for ETGs (Faber & Jackson 1976), but it has a profoundly different physical meaning. In this relation β and $L_0^{\prime }$ are free time-dependent parameters that can vary considerably from galaxy to galaxy, according to the mass assembly history and stellar evolution of each object. This relation empirically encrypts the effects of all the above-mentioned physical processes in terms of luminosity and velocity dispersion variations, parameters that can vary across time because galaxies evolve, merge, and interact.

In our previous works we tried to highlight some of the advantages offered by coupling the VT with the $L=L_0^{\prime }{\sigma }^{\beta }$ law. The first efforts were dedicated to understanding the origin of the fundamental plane (FP) of ETGs and the distributions observed in its 2D projections (D’Onofrio et al. 2017; D’Onofrio et al. 2019; D’Onofrio et al. 2020; D’Onofrio & Chiosi 2021). While discussing these problems, D’Onofrio et al. (2017) and D’Onofrio & Chiosi (2022) advanced the idea that the explanation invoked for the origin of the FP tilt (and its small scatter) should also account for the observed distributions of galaxies in all the 2D projections of the FP. The solution was found in the coupling of the VT with the time-dependent $L=L_0^{\prime }{\sigma }^{\beta }$ relation.

The key idea behind this approach is that the luminosity of galaxies is not simply related to the total stellar mass, but also to random variations caused by mergers and interactions. Accepting the $L=L_0^{\prime }{\sigma }^{\beta }$ law as an empirical descriptor of the possible changes occurring in σ and L, this implies that we can describe a galaxy with two different independent equations: the classical scalar VT based on the notion that galaxies are always very close to mechanical equilibrium, and the $L=L_0^{\prime }{\sigma }^{\beta }$ law that fully accounts for all possible processes taking place during the lifetime of a galaxy.

Our previous studies have successfully shown that a β parameter changing with time and assuming either positive and negative values, can easily explain the movements and distribution in the planes of the Sc–Rs. This approach is able to explain in a natural way the tilt of the FP, the existence of the zone of exclusion (ZoE) observed in many Sc–Rs, and the direction of motion derived from the changes in σ and L. In this work we aim to provide evidence that such an approach gives a global interpretation of the Sc–Rs observed for ETGs and that even the classical monolithic view of mass assembly is in agreement with the idea of a variable β parameter, and thus confirming the $L=L_0^{\prime }(t){\sigma }^{\beta (t)}$ law.

The paper is organized as follows. Section 2 gives a short description of the samples of galaxies (both real and simulated) used in this work. Section 3 is dedicated to the derivation of the new equations of galaxy evolution and to the different relations among the structural parameters in all FP projections. Section 4 presents a few new simple models of ETGs growing with a star formation rate (SFR) that depends on σ and accounts for the role of β. Finally, Sect. 6 provides our discussion and conclusions. In all calculations we used the parameters of the ΛCDM cosmology.

2. Samples of real and model galaxies

Observational data. The observational data used in this work are the same as used in D’Onofrio & Chiosi (2021, 2022). The data for the real galaxies are from the WINGS and Omega-WINGS databases (Fasano et al. 2006; Varela et al. 2009; Cava et al. 2009; Valentinuzzi et al. 2009; Moretti et al. 2014, 2017; D’Onofrio et al. 2014; Gullieuszik et al. 2015; Cariddi et al. 2018; Biviano et al. 2017).

The sample is not homogeneous because the spectroscopic database is only a subsample of the whole optical sample. The ETGs with available velocity dispersion σ, stellar mass M_s, and SFR, are less numerous than those extracted from the photometric database (e.g., providing R_e, I_e, n, L_V). In particular we used the following parameters: (1) the velocity dispersion σ of ∼1700 ETGs. The σ measurements come from the SDSS and NFPS databases (Bernardi et al. 2003; Smith et al. 2004) and were measured within a circular area of 3 arcsec around the center of the galaxies; (2) the luminosity, effective radius, and effective surface brightness in the V-band of several thousand ETGs derived by D’Onofrio et al. (2014) with the software GASPHOT (Pignatelli et al. 2006). The effective radius is determined from the luminosity growth curve by considering the circle that contains half the total luminosity. The effective surface intensity follows directly from the knowledge of L and R_e; (3) the distance of the galaxies derived from the redshift measured by Cava et al. (2009) and Moretti et al. (2017); (4) the stellar mass obtained by Fritz et al. (2007), only for the galaxies of the southern hemisphere.

The cross-match between the spectroscopic and optical samples provides here only 480 ETGs with available stellar mass, luminosity, velocity dispersion, Sérsic index, effective radius, and effective surface brightness. The error on these parameters is ≃20%. They are not shown in our plots because they are much lower than the observed range of variation in the structural parameters in the scaling relations and do not affect the whole distribution of ETGs.

Occasionally, we also used the catalog by Burstein et al. (1997), which contains globular clusters (GCs), dwarf galaxies (DGs) of different types, late- and early-type galaxies (LTGs and ETGs, respectively), and clusters of galaxies (GCGs). These objects are used to obtain a general idea of the Sc–Rs for systems of different sizes, but that are dynamically close to the virial condition. At times we also used the sample of Bernardi et al. (2010), but only for ETGs.

Simulated galaxies. Hydrodynamic simulations are probably the best galaxy models currently available to compare theory with observations even though several problems still bias their results. There are several suites of galaxy simulations in cosmological context, among which we recall Illustris-1 by Vogelsberger et al. (2014), Genel et al. (2014), Nelson et al. (2015), recently superseded by Illustris-TNG by Springel et al. (2018), Nelson et al. (2018), Pillepich et al. (2018a), and EAGLE by Schaye et al. (2015). We decided to adopt here Illustris-1 for two reasons. First, and most importantly, we want to be consistent with the results shown in our previous papers on this same subject that were based on the Illustris-1 models. Second, we checked that the main results of our analysis do not change passing from Illustris-1 to Illustris-TNG.

The kind of analysis carried out here is independent of the level of precision reached by models from different sources, because we are mainly interested in presenting a new method for deciphering the information encrypted in the observational data about the past history of ETGs. To this end, we extracted from the Illustris-TNG database at redshift z = 0 a sample of about one thousand model galaxies of all possible masses that are used to support the above statement.

Our data set extracted from Illustris-1 consists of several subsets of about ∼2400 galaxies each, sampled at different redshifts from to z = 0 to z = 4. A full description of these data is given in Cariddi et al. (2018) and D’Onofrio et al. (2020). In particular, we collected the effective radii, the total luminosities, the stellar masses, the velocity dispersions, the ages, and the star formation rates, together with the radii, masses, and velocity dispersions of the dark matter components.

A detailed analysis of the differences between Illustris-1 and Illustris-TNG data is made by Pillepich et al. (2018a,b), Rodriguez-Gomez et al. (2019), and Huertas-Company et al. (2019). One of major tensions between the two suites of models concerns the radii of the low mass galaxies, which are roughly M_s ≤ 5 × 10¹⁰ M_⊙, where the Illustris-TNG radii are about a factor of two smaller than those of Illustris-1, while above this value they are nearly equal (Pillepich et al. 2018a,b; Rodriguez-Gomez et al. 2019). Huertas-Company et al. (2019) compared the log(R_e)−log(M_s) plane built with the two sources mentioned above and the SDSS data of Meert et al. (2015), and find the same result (see their Fig. 11).

To better illustrate the difference in Fig. 1, we compared the data of Illustris-1 with those of Illustris-TNG-100 and the WINGS objects. The difference in the low mass range is confirmed, but the hockey stick shape of the distribution of model galaxies in the two samples is the same (see also Figs. 9 and 11 below).

Fig. 1. The size-mass and Faber–Jackson relationships. Left panel: stellar mass vs. radius relations for the Illustris-1 (open red squares) and the Illustris-TNG-100 (blue dots) samples at z = 0 and comparison of the models with the WINGS data (black dots). There are 2400 objects for the Illustris-1 sample and about 600 objects for the Illustris-TNG-100. The mean radii of Illustris-1 are smaller by about a factor of two for stellar masses lower than about 6 × 1010 M⊙, while they are nearly equal (if not slightly larger) above this limit. Right panel: L − σ plane of the same data. The symbols and color-coding are the same as in the left panel.

In addition, there is the claim that Illustris-1 simulations do not produce a realistic red sequence of galaxies due to insufficient quenching of the star formation with too few red galaxies (Snyder et al. 2015; Bottrell et al. 2017a,b; Nelson et al. 2018; Rodriguez-Gomez et al. 2019), while the Illustris-TNG simulations produce a much better red sequence (Nelson et al. 2018; Rodriguez-Gomez et al. 2019). There is also the problem of the insufficient number of red galaxies with respect to the observed population of ETGs. This is of little importance for our analysis because we do not make use of colors, but only of total luminosities.

For the internal structure of the Illustris-1 galaxies, Bottrell et al. (2017b) measured the Sersic index, the axis ratio, and the radii of these galaxies and found that too few bulge-dominated objects are produced, which is in tension with observations. In contrast, the Illustris-TNG galaxies have much better internal structural parameters (Rodriguez-Gomez et al. 2019). Fortunately, the point mass view of the Illustris-1 models we adopted ensures that our analysis is not affected too much by this problem.

Finally, the Illustris-1 data set does not give information about the morphology of the galaxies. This means that in our comparison ETGs and late-type objects are mixed in our plots. Again Fig. 11 of Huertas-Company et al. (2019) shows us that ETGs and late-type objects follow very similar trends in the Sc–Rs. The basic features of the Sc–Rs shown by ETGs are not destroyed with the addition of late-type objects.

In our work we do not make any predictions, but only qualitatively compare observations and simulations. Their behavior in the FP projections, the good match between data and simulations, and the fact that model galaxies are able to reproduce some particular features visible in the FP projections (e.g., the position of the BCGs and the existence of a ZoE) all lend support to the scenario proposed here. This makes us confident that the simulations produce galaxies with luminosity and primary structural parameters not too far from those of real galaxies.

Given the large heterogeneity of data used here, we note that the completeness of the data sample is not fundamental for the conclusions drawn in this work because we do not make any statistical analysis of the data nor do we fit any distributions. The data are only used qualitatively to show that our calculations are in agreement with the observed distributions of ETGs in the main Sc–Rs. The purpose of this paper is only that of proposing a new possible framework to analyze the evolution of ETGs.

3. Equations of galaxy evolution

The equations tracking the evolution of galaxies are based on two hypotheses. First, ETGs are always close to the virial equilibrium, a reasonable assumption since the dynamical timescale to reach this condition is on the order of the free-fall time (< 300 Myr). Second, the $L=L_0^{\prime }{\sigma }^{\beta }$ law mirrors the effects of many internal and external events affecting luminosity and velocity dispersion (i.e., mass). The two equations are

$$\begin{array}{cc}& {\sigma }^2=\frac{G}{k_v}\frac{M_{\mathrm{s}}}{R_{\mathrm{e}}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\sigma }^{\beta }=\frac{L}{L_0^{\prime }}=\frac{2\pi I_{\mathrm{e}}{R}_{\mathrm{e}}^2}{L_0^{\prime }},\hfill \end{array}$$(3)

where k_v is the non-homology parameter defined by Bertin et al. (2002). The unknown variables of this system of equations to be found are β and $L_0^{\prime }$.

Combining these two equations together, we can write

$$\begin{array}{c}\hfill log{R}_{\mathrm{e}}=alog\sigma +b{⟨\mu ⟩}_{\mathrm{e}}+c,\end{array}$$(4)

where the coefficients

$$\begin{array}{c}\hfill \begin{array}{cc}& a=(2+\beta )/3\hfill \\ \hfill & b=0.26\hfill \\ \hfill & c=-10.0432+0.333\ast (-log(G/k_v)-log(M/L)\hfill \\ \hfill & -2\ast log(2\pi )-log(L_0^{\prime }))\hfill \end{array}\end{array}$$(5)

are written in terms of β and $L_0^{\prime }$. The similarity with the FP equation is clear. This is the equation of a plane in the log(σ)−⟨μ⟩_e − log(R_e) space. The novelty is that each galaxy follows independently an equation like this. In this case, since β and $L_0^{\prime }$ are time dependent, the equation is telling us which is the instantaneous direction of motion of an object in the log(σ)−log(I_e)−log(R_e) space and in its projections.

Before showing this, let us look at the past history of the reasoning presented in this section. Starting from the same arguments and Eqs. (3) and (4), after tedious algebraic manipulations D’Onofrio & Chiosi (2022) arrived to a cubic equation in the variable β (their Eq. (10)), the coefficients of which where function of σ, I_e, R_e, M_s, and L. The cubic equation was applied to real galaxies of the WINGS list and model galaxies of the Illustris-1 catalog. In most cases three real roots were found, two of them positive and one negative. In some cases the solutions were complex and this was attributed to insufficient accuracy in the input parameters. The mutual agreement between the two sets of data (WINGS and Illustris-1) was considered to be a strong hint for self-consistency of the whole approach. This agreement was misleading because it masked first an algebraic mistake made while carrying out the lengthy analytical manipulations (i.e., a factor of 0.5 missing in front of a group of terms in logarithmic form), and second that the agreement between WINGS and Illustris-1 made via the cubic equation was in reality a circular argument as in each case the results would have been the same regardless of whether the equation was correct or not. Furthermore, attempts to incorporate the cubic equation in model galaxies did not lead to a clear understanding of the physical role and meaning played by the three $L=L_0^{\prime }{\sigma }^{\beta }$ relations associated with each time step (the factor $L_0^{\prime }$ being derived from the real luminosity by comparison). It was clear that some of the βs changed sign in the course of evolution, and also that complex solutions could occur during the lifetime of a galaxy, the low mass galaxies in particular. However, from these results the tantalizing suggestion emerged that a solution of the puzzle could be reached by changing strategy. All this led us to revise the whole problem, and thus to discovering the analytical mistake and putting the mathematical formulation on the right track. The new version of the problem is presented below. The cubic is replaced by a system of equations in the unknowns β and log $L_0^{\prime }$, which thus allows us to fully determine the $L=L_0^{\prime }{\sigma }^{\beta }$ and its evolutionary history.

Starting from Eqs. (3) and (4), after some algebra it is possible to write all the relations among the parameters of the FP projections. For the I_e − R_e plane we have

$$\begin{array}{c}\hfill I_{\mathrm{e}}=\mathrm{\Pi }{R}_{\mathrm{e}}^{\gamma },\end{array}$$(6)

where

$$\begin{array}{c}\hfill \gamma =\frac{(2/\beta )-(1/2)}{(1/2)-(1/\beta )}\end{array}$$

and Π is a factor that depends on k_v, M/L, β, and $L_0^{\prime }$ and is described by

$$\begin{array}{c}\hfill \mathrm{\Pi }={\left[{\left(\frac{2\pi }{L_0^{\prime }}\right)}^{1/\beta }{\left(\frac{L}{M_{\mathrm{s}}}\right)}^{(1/2)}{\left(\frac{k_v}{2\pi G}\right)}^{(1/2)}\right]}^{\frac{1}{1/2-1/\beta }}.\end{array}$$

For the R_e − σ plane we have

$$\begin{array}{c}\hfill R_{\mathrm{e}}=\left[\left(\frac{k_v}{G}\right)\left(\frac{L_0^{\prime }}{2\pi }\right)\left(\frac{1}{M_{\mathrm{s}}}\right)\left(\frac{1}{I_{\mathrm{e}}}\right)\right]{\sigma }^{(2+\beta )},\end{array}$$(7)

for the I_e − σ plane we have

$$\begin{array}{c}\hfill I_{\mathrm{e}}=\left[\left(\frac{G}{k_v}\right)\left(M_{\mathrm{s}}\right)\left(\frac{L_0^{\prime }}{2\pi }\right)\left(\frac{1}{R_{\mathrm{e}}^3}\right)\right]{\sigma }^{(\beta -2)},\end{array}$$(8)

and for the R_e − M_s plane we have

$$\begin{array}{c}\hfill R_{\mathrm{e}}={\left[\left(\frac{G}{k_v}\right){\left(\frac{L_0^{\prime }}{2\pi }\right)}^{2/\beta }{\left(\frac{1}{I_{\mathrm{e}}}\right)}^{2/\beta }\right]}^{\beta /(\beta +4)}{M}_{\mathrm{s}}^{\beta /(\beta +4)}.\end{array}$$(9)

It should be noted here that these equations do not represent the true physical link between two variables because their proportionality factor contains other variables as well. In other words, they do not tell us how R_e and I_e vary when σ changes. They are intermediate mathematical expressions yielding the structural parameter R_e or I_e as functions of the others. Figure 2 gives an idea of the degree of precision in reproducing the structural parameters when Eqs. (6)–(9) are used. The x-axis contains the measured parameters, while the y-axis the values calculated on the basis of our equations. The scatter in log units is in the range 0.3 − 0.4, so a factor of 2 − 2.5 uncertainty is possible and likely attributable to the ∼20% errors of the scaling parameters.

Fig. 2. Comparison between observed and calculated parameters. The data from the WINGS database are used. The black solid line gives the 1:1 relationship. The red solid line is the bi-linear least-squares fit of the distribution.

The importance of these equations is that, starting from them, we can also write the following equations (in log form):

$$\begin{array}{cc}& \beta [log(I_{\mathrm{e}})+log(G/k_v)+log(M_{\mathrm{s}}/L)+log(2\pi )+log(R_{\mathrm{e}})]\hfill \\ \hfill & +2log(L_0^{\prime })-2log(2\pi )-4log(R_{\mathrm{e}})=0,\hfill \end{array}$$(10)

$$\begin{array}{cc}& \beta log(\sigma )+log(L_0^{\prime })+2log(\sigma )+log(k_v/G)-log(M_{\mathrm{s}})\hfill \\ \hfill & -log(2\pi )-log(I_{\mathrm{e}})-log(R_{\mathrm{e}})=0,\hfill \end{array}$$(11)

and assuming

$$\begin{array}{cc}\hfill A& =log(I_{\mathrm{e}})+log(G/k_v)+log(M_{\mathrm{s}}/L)+log(2\pi )\hfill \\ \hfill & +log(R_{\mathrm{e}}),\hfill \\ \hfill B& =-2log(2\pi )-4log(R_{\mathrm{e}}),\hfill \\ \hfill A\prime & =log(\sigma ),\hfill \\ \hfill B\prime & =2log(\sigma )-log(G/k_v)-log(M_{\mathrm{s}})-log(2\pi )\hfill \\ \hfill & -log(I_{\mathrm{e}})-log(R_{\mathrm{e}}),\hfill \end{array}$$(12)

we can write the system

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A\beta +2log(L_0^{\prime })+B& =0\hfill \\ \hfill A\prime \beta +log(L_0^{\prime })+B\prime & =0\hfill \end{array}\end{array}$$(13)

with the following solutions:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \beta & =\frac{-2log(L_0^{\prime })-B}{A}\hfill \\ \hfill log(L_0^{\prime })& =\frac{A\prime B/A-B\prime }{1-2A\prime /A}·\hfill \end{array}\end{array}$$(14)

In other words, it is possible to derive the values of β and $L_0^{\prime }$ for each galaxy. This means that the knowledge of the structural parameters reveals the basic step of galaxy evolution encoded in the parameters β and $L_0^{\prime }$.

Figure 3 shows the histograms of the distributions of the β parameter derived for the galaxies of the WINGS and Illustris-1 samples. In the upper panel the β values for the Illustris-1 data are obtained from the ΔL/Δσ ratio measured on the L − σ plane. This is possible by considering the values of L and σ at two close redshift epochs (z = 0.2 and z = 0). In the lower panel we consider only the objects that are close to the virial equilibrium, those for which

$$\begin{array}{c}\hfill 2log(\sigma )=log(G/k_v)+log(M_{\mathrm{s}}/L)+log(2\pi )+log(I_{\mathrm{e}})+log(R_{\mathrm{e}})\end{array}$$(15)

Fig. 3. The distribution of the β parameter. Upper panel: histogram of the β solutions derived from Eq. (14) for the real ETGs (black line) compared with the distribution derived for the galaxies of the Illustris-1 simulation when β is calculated looking at the variation in luminosity and velocity dispersion in two close redshift epochs (z = 0.2 and z = 0). Lower panel: histogram of the β solutions derived from Eq. (14) for the real ETGs (black line) and the Illustris-1 galaxies (red line) that are close to the virial equilibrium. The solid black line gives the average value of β.

within a 20% uncertainty, and calculated β using our new analytical equations. When this condition is satisfied we obtain that 2A′/A = 1 and that β and $L_0^{\prime }$ diverge.

Notably, the values of β are both positive and negative and there is a clear deficiency of objects with β close to 0. This is true both for WINGS and Illustris-1. The average value of β is −2.44 with a rms scatter of ∼178. The positive values range from 1.05 to 1531, while the negative ones from −5.4 to −3860. The importance of Eq. (14) is that we have now an empirical thermometer of the virial condition, realized when β and $L_0^{\prime }$ diverge.

The meaning of Fig. 3 is that galaxies during their evolution can acquire either positive or negative values of β, depending on the particular events experienced (e.g., merging, stripping, star formation), and this has immediate effects on the structural parameters in the Sc–Rs, which change accordingly. Consequently, the Sc–Rs seen in their temporal framework become sources of information for the global evolution of the stellar systems.

Figure 4 shows the distribution of β as a function of the degree of virialization, expressed by the quantity 1 − 2A′/A derived. The large values of β are attained by objects very close to the virial condition. On the other hand, small β values belong to objects still away from this condition.

Fig. 4. The β parameter as proxy of virialization. The plot shows how β changes as a function of the degree of virialization (see Eq. (14)).

In closing this section we show the distribution of the β values we would obtain with the galaxy models of the Illustris-TNG-100 sample, and compare it with those of Illustris-1 and WINGS data. The three histograms are shown in Fig. 5. The difference is very small and mostly due to the smaller number of galaxies in the Illustris-TNG sample.

Fig. 5. Histogram of the β solutions derived from Eq. (14) for the real ETGs (black line) of WINGS, the Illustris-1 galaxies (red line) that are close to the virial equilibrium, and the Illustris-TNG galaxies (blue line) in the same condition. The three histograms are qualitatively similar, thus confirming that our analysis depends little on the choice between the two theoretical databases.

3.1. Trends in the FP projections

In this section we try to better explain the reasons why β can change sign during the life of a galaxy or when passing from one galaxy to another. The advantage of knowing β is much clearer when we look at the projections of the FP. Equations (6)–(8) can be further elaborated to eliminate the dependence on I_e and R_e present in their zero-points. We obtain

$$\begin{array}{cc}& I_{\mathrm{e}}={\left[\frac{G}{k_v}\frac{L_0^{\prime }}{2\pi }{M}_{\mathrm{s}}{\mathrm{\Pi }}^{3/\gamma }\right]}^{\frac{\beta -2}{1+3/\gamma }}{\sigma }^{\frac{\beta -2}{1+3/\gamma }},\hfill \end{array}$$(16)

$$\begin{array}{cc}& R_{\mathrm{e}}=\left[\frac{G}{k_v}\frac{L_0^{\prime }}{2\pi }\frac{M_{\mathrm{s}}}{\mathrm{\Pi }}\right]{\sigma }^{\frac{\beta -2}{3+\gamma }},\hfill \end{array}$$(17)

$$\begin{array}{cc}& R_{\mathrm{e}}={\left[{\left(\frac{G}{k_v}\right)}^{\beta /2}\frac{L_0^{\prime }}{2\pi }\frac{1}{\mathrm{\Pi }}\right]}^{\frac{2(\beta -2)}{{\beta }^2-6\beta +12}}{M}_s^{\frac{{\beta }^2-2\beta }{{\beta }^2-6\beta +12}}.\hfill \end{array}$$(18)

These relations now better represent the mutual dependence of the structural parameters (e.g., of I_e as a function of R_e and σ and of R_e as a function of σ) and clarify the role of β. When a galaxy moves in the L − σ plane, according to the values of β, it does the same in the other FP projections, according to the slopes reported in the last three columns of Table 1. These slopes depend on β and indicate the direction of motions (indicated by the arrows) that are visible in Figs. 6–11.

Fig. 6. Ie − Re plane. The black and red arrows show the direction of motion of galaxies in this plane for large negative and positive values of β. The black solid line gives the least-squares fit of the data, while the broken line represents the zone of exclusion.

Fig. 7. Ie − Re plane for the WINGS and Illustris-1 data. The black and red arrows show the direction of motion of the WINGS galaxies in this plane for large negative and positive values of β. The green and blue arrows are those of Illustris-1 for negative and positive values of β. In the plot only 1/10 of the Illustris-1 galaxies are used in order to distinguish the objects with different β values moving in different directions.

Fig. 8. Re − σ plane. Symbols and colors as in Fig. 6.

Fig. 9. Re − σ plane for WINGS and Illustris-1. Symbols and colors as in Fig. 7.

Fig. 10. Ie − σ plane. Symbols and colors as in Fig. 6.

Fig. 11. Ie − σ plane for WINGS and Illustris-1. Symbols and colors as in Fig. 7.

Figure 6 shows the case of the I_e − R_e plane. The black and red arrows give the direction of motion of galaxies predicted on the basis of their negative and positive values of β, respectively. Since the WINGS galaxies are well virialized, the values of β are always very large, either positive and negative. Both such slopes consistently give a direction of motion close to ∼ − 1 in the I_e − R_e plane. The −1 slope is that predicted on the basis of the VT (represented by the broken line, which is also the ZoE; see D’Onofrio & Chiosi 2021). We note that no galaxies can cross the ZoE because their motion is nearly parallel to the ZoE.

These arguments demonstrate the reason why there is a ZoE in the I_e − R_e plane: the only possible direction of motion for well virialized objects is that with slope ∼ − 1 generated by the large positive and negative values of β. When β ∼ 0 (i.e., when the galaxies are less virialized), they can move in other directions in this plane. Unfortunately, our sample does not include the dwarf ETGs, which are usually distributed in a cloud, below the ZoE with radii lower than 3−4 kpc (see, e.g., Capaccioli et al. 1992; D’Onofrio et al. 2020). For these objects we predict values of the slopes in all possible directions.

The only way to check this is to make use of the model galaxies either of Illustris-1 or Illustris-TNG. Figure 7 confirms our prediction. Although the data of Illustris-1 are affected by the well-known problem of systematically larger R_e with respect to the observed values (D’Onofrio et al. 2020; Bottrell et al. 2017b), we note that several objects have arrows nearly orthogonal to those of the well-virialized galaxies. The expected motions of the dwarf galaxies are in all possible directions, thus giving rise to the cloud of the “ordinary” ETGs defined by Capaccioli et al. (1992). Furthermore, when the curve of the distribution changes, we note a progressive variation in the arrow directions. This means that the overall distribution in this plane is governed by the different movements of the galaxies in the L − σ plane described empirically by the different values of β and $L_0^{\prime }$.

The direction of the arrows displayed in each figure visualizes the expected displacement of a galaxy based on the actual value of β. However, the arrows only give the direction of motion, not the orientation of the future temporal evolution of a galaxy. Furthermore, they do not indicate the path followed by each galaxy to reach the current observed position in the diagrams.

The same can be said for the other two FP projections. Figures 8 and 10 represent the R_e − σ and I_e − σ planes, respectively. Here the role of β is much clearer. It is evident that the galaxies with negative β values move in different directions with respect to those with positive β, originating the curvatures observed in these diagrams.

Once more the addition of the Illustris-1 data (Figs. 7 and 9) confirms that the slopes derived from the βs are consistent with the observed distribution of ETGs and demonstrates that the observed curvatures originate from the different motion of galaxies with positive and negative values of β.

Figures 12 and 13 are similar plots for the log(R_e)−log(M_s) plane. Even in this important diagram we observe a ZoE (dashed line with slope equal to 1). Our calculations predict why this ZoE is here: all the virialized objects (with large β values) can only move in the direction with slope equal to 1 (see Table 1). The R_e of Illustris-1 are somewhat larger than those measured, but the general behavior is in good agreement with the observed distribution. The galaxies with large positive and negative values of β move with a slope close to 1, while in the cloud of points with low masses we can find objects with different directions.

Fig. 12. The log(Re)−log(Ms) plane for the WINGS galaxies. Symbols and colors as in Fig. 6.

Fig. 13. The log(Re)−log(Ms) plane for WINGS and Illustris-1 galaxies. Symbols and colors as in Fig. 7.

We conclude that all the positions of ETGs in the FP projections and in the log(R_e)−log(M_s) plane depend on the motions that occurred during the peculiar evolutionary path followed by each galaxy. When the galaxies are well virialized these motions can occur only in fixed directions depending on the value of β. This is a coherent and self-consistent explanation of all the main scaling relations built with the structural parameters of ETGs. It follows that even the FP, the “father” of the scaling relations, should be explained in a similar way.

Looking at Table 1 in detail we also note that,

– in all FP projections, when β becomes progressively negative (i.e., when the objects rapidly decrease in luminosity at nearly constant σ), the slopes either converge to the values predicted by the VT (in the I_e − R_e relation and in the R_e − M_s relation) or diverge toward large values (in the I_e − σ and R_e − σ relations) because the galaxy keeps its velocity dispersion when the luminosity decreases (only I_e and R_e vary). This offers a natural explanation of the ZoE.

– The table also suggests that positive and negative values of β are equally permitted with both real and simulated data. In general, the objects that are still active in their star formation or have recently experienced a merger have positive values of β, while those progressively quenching their SF have increasing negative β.

– The curvature in the observed distributions (i.e., the transition from the large cloud of small galaxies to the much narrower tail of the brightest objects) is naturally explained by the existence of positive and negative values of β.

A way to better understand the effects played by β is to think of the possible variations in R_e and I_e when L and σ vary in the L − σ plane. There are four possible changes of L and σ in this plane. They are schematically shown in Table 2, which displays, according to the values of β, the expected variations in R_e and I_e, when L, M_s, and σ vary. We note that when β is negative there is not necessarily a decrease in luminosity, and when β is positive a decrease in luminosity might also occur.

When the luminosity of a galaxy changes, both the effective radius and the mean effective surface intensity I_e vary. This happens because R_e is the radius of the circle that encloses half of the galaxy total luminosity. Since the ETGs have different stellar populations with different ages and metallicity, it is highly improbable that the decrease in luminosity does not change the whole appearance of the luminosity profile². Consequently the growth curve changes and determines a variation in R_e and I_e. If the luminosity decreases passively, in general we could expect a decrease in R_e and an increase in I_e. On the other hand, if a shock induced by harassment or stripping induces an increase in L (and a small decrease in σ), we might expect an increase in R_e and a decrease in I_e.

The observed variations of these parameters depend strongly on the type of event that a galaxy is experiencing (stripping, shocks, feedback, merging, etc.). In general, one should keep in mind that these three variables L, R_e and I_e are strongly coupled each other and that even a small variation in L might result in ample changes in R_e and I_e. In this context we begin to understand that the Sc–Rs are useful tools for guessing both the dynamics and the evolutionary state of the stellar content of a galaxy.

In summary, what we claim here is that all the above diagrams should be analyzed taking into account the effects of time, and should not be investigated separately. They are snapshots of an evolving situation, and this temporal evolution cannot be discarded. The $L=L_0^{\prime }{\sigma }^{\beta }$ law catches such evolution in the correct way by predicting the direction of the future motion of each galaxy in the diagnostic planes (D’Onofrio & Chiosi 2023). In principle, this way of reasoning should allow us to understand why galaxies are in the positions observed today in each diagram. As β gives only the present direction of motion and not that of the motion in the past, the simultaneous use of simulations and high redshift observations might help us to infer the possible precursors of the present-day galaxies on the basis of the physical properties and the distribution in the FP projections, indicated by the values of β. In other words, these scaling relations become a possible tool for inferring the evolutionary path of each galaxy.

3.2. Origin of the FP and its tilt

The final step is related to the question of the origin of the FP. Equation (4) tells us that each galaxy follows its own FP-like equation, whose coefficients are functions of β. Starting from Eq. (4), it is possible to derive the coefficients a, b, and c of the plane hosting each single ETG. To do this we adopt the notation that is commonly used for the FP, in which ⟨μ⟩_e is expressed in mag arcsec⁻² and R_e in kpc:

$$\begin{array}{c}\hfill log(R_{\mathrm{e}})=alog(\sigma )+b{⟨\mu ⟩}_{\mathrm{e}}+c,\end{array}$$(19)

where a, b, and c are from Eq. (4). This transformation is necessary because in our notation I_e is expressed in L_⊙ pc⁻² and R_e is in pc.

The distribution of these coefficients for all the WINGS sample of 479 galaxies is shown in Fig. 14. It is clear from the figure that the values for the FP coefficients, derived from the fit of the ETGs distribution and indicated by the dashed areas in each panel, are very close to the average of the single coefficients calculated by Fig. 14 (vertical black lines). The gray bands show the interval of a, b, and c obtained by D’Onofrio et al. (2008) fitting the FP of ETGS separately for each cluster of the WINGS data set.

Fig. 14. Histograms of the values of the calculated FP coefficients a (top), b (middle), and c (bottom). In each panel the histogram of the coefficient (blue line), the average value (vertical black line), and the median (vertical red line) are shown. In the case of the b coefficient, which does not depend on β, the average and the median are the same (the black and red vertical lines coincide). The dashed gray regions give the intervals of the FP coefficients found by fitting the distribution of ETGs in the log(σ)−log(Ie)−log(Re) space.

In other words, we have framed the FP and its tilt in a new context where each ETG follows its own Eq. (4), namely FP, and contributes to shaping the global FP (both tilt and thickness) of the ETG population. Since the FP coefficients are obtained from a fit, it is clear that the final coefficients of the plane will be close to the average of the single values valid for each object. Some differences are expected because the final values depend on the sample adopted (each having its own average) and on the technique used to perform the fit.

With this statement we do not mean to say that the various mechanisms invoked in the literature to explain the tilt and thickness are incorrect. Rather, we claim that all of them can actually contribute to the average properties of the galaxy sample, giving rise to a different β for each object.

Since its discovery, the FP has been the subject of several studies aimed at understanding why the plane is tilted with respect to the prediction of the VT (M_s ∝ R_eσ²), and why its intrinsic scatter is so small (see, e.g., Faber et al. 1987; Ciotti 1991; Jorgensen et al. 1996; Cappellari et al. 2006; D’Onofrio et al. 2006; Bolton et al. 2007, among many others). While the VT predicts a = 2 and b = −1, the values coming from the fit of several samples of ETGs are systematically lower (a ∼ 1.2) and higher (b ∼ −0.8), and vary according to the sample and the fitting strategy used.

Among the physical mechanisms invoked to explain the FP tilt we find the following: (1) a progressive change in the stellar mass-to-light ratio (M_s/L), (see, e.g., Faber et al. 1987; van Dokkum & Franx 1996; Cappellari et al. 2006; van Dokkum & van der Marel 2007; Holden et al. 2010; de Graaff et al. 2021); (2) structural and dynamical nonhomology (see, e.g., Prugniel & Simien 1997; Busarello et al. 1998; Trujillo et al. 2004; D’Onofrio et al. 2008); (3) dark matter (DM) content and distribution (see, e.g., Ciotti et al. 1996; Borriello et al. 2003; Tortora et al. 2009; Taranu et al. 2015; de Graaff et al. 2021); (4) star formation history (SFH) and initial mass function (IMF; see, e.g., Renzini & Ciotti 1993; Chiosi et al. 1998; Chiosi & Carraro 2002; Allanson et al. 2009); (5) the effects of environment (see, e.g., Lucey et al. 1991; de Carvalho & Djorgovski 1992; Bernardi et al. 2003; D’Onofrio et al. 2008; La Barbera et al. 2010; Ibarra-Medel & López-Cruz 2011; Samir et al. 2016). In addition, recent observational work has shown that variations in the M_s/L ratio can account for only half of the tilt (see D’Eugenio et al. 2021), with the remainder being due to structural variation and possibly variations in the galaxy-averaged initial mass function of the stellar populations. Uncertainties in M_s can affect the tilt if the error is mass-dependent, although this systematic uncertainty is not large enough (see Leja et al. 2019; Lower et al. 2020). Schechter et al. (2014), using strong lensing measurements, provide an independent estimate of M_s but still find a tilt of the FP. Finally, the tilt is found in cosmological simulations (Rosito et al. 2019a,b; de Graaff et al. 2023).

All these effects are in practice involved in our view of the problem. Since each sample of galaxies has its own average value of β, because of the different history of mass accretion and luminosity evolution, it is easy to verify that systematic changes in the tilt could arise for the above-mentioned reasons. When the sample changes its average properties, a small variation in the tilt of the FP follows. This explains for instance why Robertson et al. (2006) found that star-forming and quiescent galaxies follow different Sc–Rs (i.e., a different FP tilt), due to differences in the merger histories.

With the emerging of the hierarchical scenario of galaxy formation and evolution, some additional mechanisms for the FP tilt have been proposed: (1) the effects of dissipationless merging (Nipoti et al. 2003); (2) gas dissipation (Robertson et al. 2006); (3) the irregular sequence of mergers with progressively decreasing mass ratios (Novak 2008); (4) multiple dry mergers of spiral galaxies (Taranu et al. 2015). Since the galaxy properties change with time, the slope of the FP is expected to change with redshift. This is confirmed by the numerical models of single galaxies, large-scale cosmological simulations, and observational surveys at different redshifts, among others (see Beifiori et al. 2017; Rosito et al. 2019a,b; Lu et al. 2019; Ferrero et al. 2021; de Graaff et al. 2023 and references therein).

The most remarkable physical feature of the FP is the observed very small scatter, which amounts to ≈0.05 dex in the V-band. It seems to require a sort of fine tuning among different physical processes. The scatter has been attributed to (1) the variation in the formation epoch; (2) the DM content; (3) the existence of metallicity or age trends; and (4) the variations in the mass-to-light ratio M/L (see, e.g., Faber et al. 1987; Gregg 1992; Guzman et al. 1993; Forbes et al. 1998; Bernardi et al. 2003; Reda et al. 2005; Cappellari et al. 2006; Bolton et al. 2008; Graves et al. 2009; Graves & Faber 2010; Auger et al. 2010; Magoulas et al. 2012). Our approach cannot predict the scatter around the FP, because this does not depend on the structural parameters, but on the properties of the stellar populations and the peculiar history of mass accretion or stripping.

Although investigating the causes of the tilt and the small dispersion around the FP is beyond the aims of this paper, we conclude this section with one consideration. We note that going back to high redshift, the numerical simulations of Illustris-1 and Illustris-TNG show that the FP and the tail of the MRR persist until redshift z ≃ 1.6 and then disappear or are no longer as well defined (Chiosi et al., in prep.). The concomitant appearance of the MRR and FP for ETGs more massive than about 10¹⁰ M_⊙ may be a mere coincidence, but it is a question that must be investigated.

The galaxies on the tails of the FP and the MRR are massive ETGs, whose stellar content is predominantly made up of old stars; in the case of mergers with objects of lower mass, the percentage of younger stars does not significantly alter the luminosity and the colors of the basic stellar populations. To quantify this statement we make use of the following example. At proceeding galaxy building via the hierarchical scenario, the probability that a massive object merges with another of similar mass becomes rarer and rarer as galaxies become more massive. Therefore, massive ETGs tend to evolve in isolation or by merging with objects of much lower mass. In general, the merger of two galaxies with very different masses (e.g., M₁/M₂ ≃ 1/10) and some companion stellar activity leaves the mass and velocity dispersion nearly unchanged, while the luminosity first undergoes a burst of short duration and relative intensity proportional to the luminosity ratio L₁/L₂ (only slightly higher than the previous value). This should correspond to a nearly vertical shift on the FP of small amplitude. The M/L ratio either remains unchanged or slightly decreases, thus causing a little scatter of the FP.

The opposite should occur in a merger between two galaxies of nearly equal mass, a typical situation in the range of low mass galaxies. In this case the mass and luminosity both change. If additional star formation occurs, there should be an additional increase in the total luminosity that depends on the amount of mass converted into new stars. Therefore, the total luminosity should hardly recover the pre-burst value, the mass-weighted mean of the two component galaxies (see Fig. 11 in Tantalo & Chiosi 2004a). So most likely the luminosity remains higher than before, and the M/L ratio is expected to decrease. This should generate a tilt of the FP in the correct direction. It is not easy to foresee the effect on the scatter. Better estimates require numerical simulations of bursts of star formation. In any case, after a short time interval the maximum shift in luminosity cannot overcome a factor of ∼2 (0.3 dex).

4. Application to model galaxies

To lend support to the picture outlined above concerning the physical role and meaning of the $L=L_0^{\prime }{\sigma }^{\beta }$ relation and the role of the parameters β and $L_0^{\prime }$ without resorting to the numerical simulations of Illustri-1 and Illustris-TNG over which we have no control at all, we make use of very simple, almost analytical models of galaxy formation and evolution. The ideal models of this type suited to describe ETGs are those in which the total mass increases by infall and the stars are formed according to a simple law of SFR developed long ago by Chiosi (1980) and extended by Tantalo et al. (1998a). The novelty here is that we incorporated the equations for β and $L_0^{\prime }$ (Eq. (4)) into the models once the luminosity, the radius, and the velocity dispersion were calculated. With the aid of these models and Eqs. (16) and (17) we calculated the basic relationships I_e − σ and R_e − σ, and finally made a cross-test of mutual consistency between the results from the galaxy models and the β and $L_0^{\prime }$ theory. These simple models of galaxy formation and evolution were first proposed by Chiosi (1980), much later extended by Tantalo et al. (1998a), and recently used by Chiosi et al. (2017) to study the cosmic SFR and by Sciarratta et al. (2019) to investigate the galaxy color-magnitude diagram. Although they may look too simplistic compared to the numerical models of Illustris-1 and Illustris-TNG, they do catch the main features of these galaxies and are suitable for our purposes.

In brief, a galaxy of total mass M_G is made of baryonic (B) and dark matter (D), of mass M_B and M_D, respectively, and at any time satisfies the following equation:

$$\begin{array}{c}\hfill M_{\mathrm{G}}(t)=M_{\mathrm{B}}(t)+M_{\mathrm{D}}(t).\end{array}$$(20)

At all times M_B(t) and M_D(t) are in cosmological proportion, which means they satisfy the condition M_D(t) = f_cM_B(t), where f_c depends on the adopted ΛCDM cosmological model of the Universe (f_c ≃ 6.1 in our case).

The baryonic mass is assumed to be originally in the form of gas, to flow in at a suitable rate, and to transform into stars when the physical conditions allow it. At the same rate dark matter is also allowed to flow in together with the baryonic matter to build up the total gravitational potential. Suitable prescriptions of their spatial distribution are needed to calculate the gravitational potential (see Tantalo et al. 1998a, for more details).

This kind of galaxy model is called an infall model, the essence of which resides in the gas accretion into the central region of the proto-galaxy at a suitable rate (driven by the timescale τ) and in the gas consumption by a Schmidt-like law of star formation. The gas accretion and consumption coupled together provide a time dependent SFR closely resembling that resulting from N-body simulations (e.g., Chiosi & Carraro 2002; Merlin & Chiosi 2006, 2007; Merlin et al. 2012).

At any time t the baryonic mass M_B is given by the sum

$$\begin{array}{c}\hfill M_{\mathrm{B}}(t)=M_{\mathrm{g}}(t)+M_{\mathrm{s}}(t),\end{array}$$(21)

where M_g(t) is the gaseous mass and M_s(t) the mass in stars. At the beginning, both the gas and the star mass in the proto-galaxy are zero M_g(t = 0) = M_s(t = 0) = 0. The rate of baryonic mass accretion (and in turn gas accretion) is driven by the timescale τ according to

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_{\mathrm{B}}(t)}{\mathrm{d}t}=M_{\mathrm{B},\tau }exp(-t/\tau ),\end{array}$$(22)

where M_B, τ is a constant with the dimensions of [Mass/Time] to be determined by imposing that at the galaxy age T_G the total baryonic mass of the galaxy M_B(T_G) is reached:

$$\begin{array}{c}\hfill M_{\mathrm{B},\tau }=\frac{M_{\mathrm{B}}(T_{\mathrm{G}})}{\tau [1-exp(-T_{\mathrm{G}}/\tau )]}·\end{array}$$(23)

Therefore, by integrating the accretion law, the time dependence of M_B(t) is

$$\begin{array}{c}\hfill M_{\mathrm{B}}(t)=\frac{M_{\mathrm{B}}(T_{\mathrm{G}})}{[1-exp(-T_{\mathrm{G}}/\tau )]}[1-exp(-t/\tau )].\end{array}$$(24)

Since dark matter flows in at the same rate of the baryonic matter, it obeys similar equations in which M_B is replaced by M_D. However, since at any time M_B and M_D, are in cosmic proportions, M_D = f_cM_B, the equations for M_D are superfluous and the normalization on M_B is enough. The underlying hypothesis is that the presence of dark matter does not affect the evolution of the baryonic component, but for its effect on the gravitational potential energy. To this end, some assumptions about the spatial distribution of M_B and M_D are needed. In other words, assuming spherical symmetry, the radii R_B and R_D must be specified.

The timescale τ is related to the collapse time and the average cooling rate of the gas. Therefore, it is expected to depend on the mass of the system. At the same time, the gas mass increases by infall and decreases by star formation.

The rate of star formation is modeled throughout the whole life of the galaxy with the Schmidt (1959) law:

$$\begin{array}{c}\hfill \mathrm{\Psi }(t)\equiv \frac{\mathrm{d}{M}_{\mathrm{s}}}{\mathrm{d}t}=\nu M_{\mathrm{g}}{(t)}^k,\end{array}$$(25)

where k regulates the dependency of the SFR on the gas content; we assume k = 1. The quantity ν is the efficiency parameter of the star formation process that must be specified (see below).

In the infall model, because of the interplay between gas accretion and consumption, the SFR starts low, reaches a peak after a time approximately equal to τ, and then declines. The functional form that could mimic this behavior is the time-delayed exponentially declining law:

$$\begin{array}{c}\hfill \mathrm{\Psi }(t)\propto \frac{t}{\tau }exp(-\frac{t}{\tau })·\end{array}$$(26)

The Schmidt law in Eq. (25) is therefore the link between gas accretion by infall and gas consumption by star formation.

As a whole, this kind of approach stands on a number of observational and theoretical arguments among which we recall first that the parameters ν and τ can be related to morphology (Buzzoni 2002) and to the presence of ongoing star formation activity inside observed galaxies (Cassarà et al. 2016), and second that these quantities can be easily tuned in order to fit observational data. Complex phenomena that would affect the rate of gas cooling, such as active galactic nuclei (AGN), can also be empirically taken into account without going into detail (see, e.g., Chiosi et al. 2017).

The infall models we describe may include many important physical phenomena, for instance gas heating by supernova explosions (both type Ia and type II), stellar winds, gas cooling by radiative emission, and the presence of galactic winds (see Tantalo et al. 1998a for details).

4.1. Outline of the galaxy models

The complexity of real globular clusters, galaxies, and galaxy clusters and the history of their evolution are reduced here to ideal systems of which we know the current masses M(t), M_B(t), M_D(t), M_s(t), and M_g(t), together with the mass abundances of some important elements X_i(t), (where i stands for H, He, C, N, O, Mg, ...Fe) and total abundance of heavy elements Z(t)³, and finally half-stellar mass (half-light) radius R_e(t) and dark mass radius R_D(t). At each time the system contains many stellar populations of different metallicities and ages that can be approximated by single stellar populations (SSPs) of mean metallicity ⟨Z(t)⟩ and mean age T(t) defined by the relation T(t) = M_s(t)/⟨Ψ(t)⟩, where ⟨Ψ(t)⟩ is the mean SFR in the interval 0 ÷ t (with t the current age). This value of the age T(t) is used to infer the current luminosity associated with the stellar content M_s(t), (see below).

The infall model of a galaxy must be completed with the radii R_e(t) and R_D(t), (which are necessary to calculate the velocity dispersion of the stellar component) and the gravitational potential for the onset of galactic winds. To this end, we briefly discuss a few items of interest.

(i) The M_D − M_L and R_D − R_L relationships. Following Bertin et al. (1992) and Saglia et al. (1992), the spatial distribution of the dark component with respect to the luminous component in dynamical models is such that the mass and radius of the dark component (M_D, R_D) are related to those of the luminous component (M_L, R_L) by

$$\begin{array}{c}\hfill \frac{M_{\mathrm{L}}(t)(t)}{M_{\mathrm{D}}(t)}\ge \frac{1}{2\pi }\frac{R_{\mathrm{L}}(t)}{R_{\mathrm{D}}(t)}[1+1.37\frac{R_{\mathrm{L}}(t)}{R_{\mathrm{D}}(t)}],\end{array}$$(27)

where we can pose M_L(t)≃M_s(t), M_D(t) = f_cM_B(t)≥f_cM_s(t) and R_L(t)≃2R_e(t). Therefore, knowing M_D(t), M_s(t), and R_e(t), we can obtain an estimate of R_D(t) to be used in the calculation of the total gravitational potential. According Bertin et al. (1992) and Saglia et al. (1992), typical values are M_L/M_D ≃ 0.2 and R_L/R_D ≃ 0.2. Consequently within R_e the mass of dark matter is low with respect to the stellar mass and can be neglected. Furthermore, the binding gravitational energy of the gas and stars is given by

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_j(t)=-{\alpha }_{\mathrm{L}}G\frac{M_j(t)M_{\mathrm{L}}(t)}{R_{\mathrm{L}}(t)}-G\frac{M_j(t)M_{\mathrm{D}}}{R_{\mathrm{L}}(t)}\mathrm{\Omega }{\prime }_{\mathrm{LD}},\end{array}$$(28)

where j stands for g (gas) or s (stars), and α_L is a numerical factor = 0.5, and finally the term

$$\begin{array}{c}\hfill \mathrm{\Omega }{\prime }_{\mathrm{LD}}=\frac{1}{2\pi }\left(\frac{R_{\mathrm{L}}(t)}{R_{\mathrm{D}}}\right)[1+1.37\left(\frac{R_{\mathrm{L}}(t)}{R_{\mathrm{D}}}\right)]\end{array}$$(29)

is the contribution to the gravitational energy given by the presence of dark matter. With the assumed ratios M_L/M_D and the above replacements of M_L, R_L, and R_D, the term ${\mathrm{\Omega }}_{\text{LD}}^{\prime }$ is about 0.04. Therefore, in the evaluation of the velocity dispersion of the stellar component via the VT the effect of DM can be neglected.

(ii) Velocity dispersion. The velocity dispersion of an object with M_D(t), M_s(t), and radius R_e(t) is derived from the scalar VT. At each time an object is assumed to be very close to the condition of mechanical equilibrium, and hence to satisfy the relation

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{s}}(t)=\sqrt{\frac{G}{k_v}\frac{M_{\mathrm{s}}(t)}{R_{\mathrm{e}}(t)}}·\end{array}$$(30)

(iii) The R_e(M_s) relation. The mass–radius relation (MRR) suited to our models is the empirical law proposed by Fan et al. (2010) in the context of the ΛCDM cosmology. The expression is

$$\begin{array}{c}\hfill R_{\mathrm{e}}=0.9\left(\frac{S_{\mathrm{S}}(n_{\mathrm{S}})}{0.34}\right)\left(\frac{25}{m}\right){\left(\frac{1.5}{f_{\sigma }}\right)}^2{\left(\frac{M_{\mathrm{D}}}{{10}^{12}{M}_{\odot }}\right)}^{1/3}\frac{4}{(1+z_f)},\end{array}$$(31)

where M_D, M_s, and R_e have their usual meaning; R_e is in kpc; z_f is the redshift at which the collapse took place; S_S(n_S) indicates the shape of the baryonic component, which in turn is related to the Sérsic brightness profile from which R_e is derived; n_S is the Sérsic index; f_σ is the three-dimensional stellar velocity dispersion as a function of the DM velocity dispersion, σ_s = f_σσ_D; and m is the ratio M_D/M_s. We adopt here S_S(n_S) = 0.34 and f_σ = 1 (for more details, see Fan et al. 2010; Chiosi et al. 2020, and references therein). The most important parameter of Eq. (31) is the ratio m = M_D/M_s, which is briefly discussed below.

The MRR of Eq. (31) is the locus of galaxy models on the MR-plane, the formation of which occurred at redshift z_f. It represents the position of model galaxies for different sources (Chiosi & Carraro 2002; Merlin et al. 2012; Vogelsberger et al. 2014); however, it does not correspond to the real MRR observed for objects from GCs to ETGs and GCGs because cosmological effects are also present (the subject is thoroughly discussed by Chiosi et al. 2020, to whom we refer for all details).

(iv) The M_D/M_s ratio. Based on the Illustris-1 data Chiosi et al. (2020) investigate how this ratio varies in the mass interval 10^8.5 < M_D < 10^13.5 (masses are in M_⊙) and from z = 0 to z = 4, and proposed the following relation:

$$\begin{array}{c}\hfill m\equiv \frac{M_{\mathrm{D}}}{M_{\mathrm{s}}}=(-0.223z_f+0.375)log{M}_{\mathrm{D}}+(3.138z_f-3.430).\end{array}$$(32)

In the present study, however, we follow a different strategy that at each time step tightly correlates the mass in stars M_s(t) to the total baryonic mass M_B(t) and to the total mass of dark matter M_D(t). At each time we have M_D(t) = f_c M_B(t), where f_c is the cosmic proportion (f_c ≃ 6). The mass in stars M_s(t) is determined by the efficiency of star formation and in any case it is a fraction of the current baryonic mass. Therefore, the parameter m is given by the relation

$$\begin{array}{c}\hfill logm(t)=log\frac{M_{\mathrm{D}}(t)}{M_{\mathrm{s}}(t)}·\end{array}$$(33)

At the beginning of a galaxy history the ratio m is very high and then declines, tending to the limit value f_c ≃ 6 if the total baryonic mass is eventually turned into stars. Examples of the time behavior of the ratio m are shown below when we present our model galaxies in detail.

(v) The star formation rate. Thanks to the short timescale of the energy input from massive stars (a few million years), compared to the mass accretion timescale by infall (from hundreds to thousands of million years) the galaxy is not supposed to differ from an equilibrium state so that the Talbot & Arnett (1975) formalism can be applied. Chiosi (1980) and Chiosi & Matteucci (1980) adopted the SFR of Talbot & Arnett (1975) to model disk galaxies in which the surface mass density of stars, gas, and total baryonic mass are used and a suitable radial distance $\stackrel{\sim }{r}$ is introduced.

We adapted their formalism to our case (in which spherical symmetry is implicitly assumed),

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_{\mathrm{s}}(r,t)}{\mathrm{d}t}=-\frac{\mathrm{d}{M}_{\mathrm{g}}(r,t)}{\mathrm{d}t}=\stackrel{\sim }{\nu }{\left[\frac{M(r,t)M_{\mathrm{g}}(r,t)}{M(\stackrel{\sim }{r},t)}\right]}^{\kappa -1}{M}_{\mathrm{g}}(r,t),\end{array}$$(34)

where M_g(r, t) and M_s(r, t) are respectively the mean mass density of gas and stars within the generic sphere of radius r at time t, $M(\stackrel{\sim }{r},t)$ is the total mass density (gas and stars) within a particular radial distance from the galaxy center, and $\stackrel{\sim }{\nu }$ is a parameter measuring the efficiency of star formation. The radius $\stackrel{\sim }{r}$ is a suitable radial scale controlling star formation. In Larson’s view they might be associated with the radial distance at which the central spheroidal component exerts its tidal effect on the residual external gas. As a consequence, at any time the SFR is significantly inhibited at distances $r>\stackrel{\sim }{r}$. Since our models do not include any geometrical description, but deal with a galaxy as a point-mass entity whose mass varies with time, we drop the radial dependence of the SFR and the rate of star formation is simply reduced to

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_{\mathrm{s}}(t)}{\mathrm{d}t}=-\frac{\mathrm{d}{M}_{\mathrm{g}}(t)}{\mathrm{d}t}=\stackrel{\sim }{\nu }{\left[\frac{M(t)M_{\mathrm{g}}(t)}{M(t)}\right]}^{\kappa -1}{M}_{\mathrm{g}}(t)=\stackrel{\sim }{\nu }{M}_{\mathrm{g}}{(t)}^{\kappa }.\end{array}$$(35)

Since in all infall models M_g(t) increases by infall and decreases by star formation, the SFR starts low, reaches a peak after a time approximately equal to τ, and then declines. By varying τ (timescale of the galaxy formation process) we can recover all types of star formation indicated by observational data, from GCs to LTGs and ETGs. The infall scheme and companion SFR are widely used in many studies on the subject of galactic chemical evolution (e.g., Matteucci 2016, for a review and references). The infall galaxy model is very flexible and can be adapted to a wide range of astrophysical problems. We recall that it was adopted by Bressan et al. (1994) to model the spectro-photometric evolution of ETGs reduced to point mass objects, extended by Tantalo et al. (1998b) to the case of spherical systems made of BM and DM mimicking ETGs, adapted by Portinari & Chiosi (2000) to include radial flows of gas in disk galaxies, and recently used by Chiosi et al. (2017) to study the cosmic SFR and by Sciarratta et al. (2019) to investigate the color-magnitude diagram of galaxies in general.

(vi) The SF efficiency $\stackrel{\sim }{\nu }$. In most galaxy models of this kind the specific efficiency of star formation $\stackrel{\sim }{\nu }$ is an external free parameter to be adjusted according to the case under investigation. In this paper we follow a different strategy and derive $\stackrel{\sim }{\nu }$ from other properties of the models. Starting from the idea put forward by Brosche (1970, 1973) that the efficiency of star formation is driven by the velocity dispersion, we suppose that $\stackrel{\sim }{\nu }$ can be written as

$$\begin{array}{c}\hfill \stackrel{\sim }{\nu }={\nu }_0{[\frac{{\sigma }_{\mathrm{t}}}{{\sigma }_{\mathrm{s}}}\times \frac{{\sigma }_T}{{\sigma }_{\mathrm{s}}}]}^{0.5},\end{array}$$(36)

where σ_s, σ_t, and σ_T are the velocity dispersion calculated using only the stellar component and the total mass, both measured at current time t and present-day age T. The factor ν₀ depends on the choice made for κ, and ensures the correct dimensions to $\stackrel{\sim }{\nu }$. For κ = 1, $\stackrel{\sim }{\nu }\equiv 1/t$. Finally, the harmonic mean between two different normalizations is meant to deal with the uncertainty affecting the whole procedure. Since the stellar mass M_s and radius R_e grow with time, the efficiency is high at young ages and decreases with time toward the limit value of ν ≃ 1.

(vii) Luminosity and specific intensity from mean SSPs. In order to calculate the B and V luminosities and the associated specific intensities I_eB and I_eV of the stellar content of galaxy models in the course of their evolution, we make use of the SSPs with the Salpeter (1955) IMF: slope in number x = −2.35, lower mass M_l = 0.1 M_⊙, upper mass M_u = 100 M_⊙, total SSP mass M_SSP = 5.82 M_⊙, metallicity from Z = 0.0004 to Z = 0.04 (six values in total), and age from 10 Myr to 14 Gyr of the library by Bertelli et al. (2008, 2009), Tantalo (2005). The absolute M_B and M_V magnitudes can be plotted against the logarithm of the age in years, and for each passband the mean age–magnitude relation is derived. Owing to the nearly linear behavior of each relationship, a linear fit is suited to obtain the relation between the mean absolute magnitude and the age t. These are given by

$$\begin{array}{cc}& M_{\mathrm{B}}=2.361logt-17.841,\hfill \end{array}$$(37)

$$\begin{array}{cc}& M_V=1.975logt-14.886.\hfill \end{array}$$(38)

The age is expressed in years. The B and V magnitudes of the original SSPs with different metallicity are shown in Fig. 15 together with the metallicity averaged SSP (full dots). The mean values of the magnitudes are meant to mimic the mixture of chemical compositions in a galaxy. At each time we know the total mass made of stars of different ages and chemical compositions. In practice, we assume that this complex situation can be reduced to a single SSP of the same mass, mean chemical composition (metallicity), and mean age T. The mean age is evaluated from the relation T(t) = M_s(t)/⟨Ψ(t)⟩, where ⟨Ψ(t)⟩ is the mean SFR in the interval 0 ÷ t (with t the current age). Using the mean age T(t), from Eqs. (37) and (38) we derive the B/V magnitudes (the luminosities) per unit mass of the SSP and then re-scale them to the mass M_s of the galaxy.

Fig. 15. Magnitude vs. age relationships for SSPs of different metallicity (in different colors): MB (top) and MV (bottom). From top to bottom the metallicity is Z = 0.0001, 0.001, 0.010, 0.019, 0.040, and 0.070. The black dotted lines are the mean values of MB and MV over the metallicity.

(viii) Solution of the basic equation Eq. (14). At each time step of the evolutionary history of a model galaxy, when the star mass M_s, the radius R_e, the velocity dispersion σ_s, the luminosities L_B and L_V (in solar units), and the specific intensities I_eB and I_eV are known, the equation system Eq. (14) is solved deriving β and $L_0^{\prime }$ at each time step. These are the two physical quantities that in our view drive the distribution of galaxies in the space of the physical parameters L, R_e, σ, and I_e, and determine the observed FP.

Final remarks. The model age refers to the galaxy rest-frame and goes from T_g = 0 at redshift z_f (when the galaxy is assumed to form) to T_g = T_G at z = 0 (present time). The corresponding ages of the Universe T_U(z) are T_U(z_f) and T_U(0). For the ΛCDM cosmology with H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.71, Ω_m = 0.23, Ω_Λ = 0.73, Ω_mD/Ω_mB ≃ 6, we obtain T_U(z_f) = 0.484 Gyr for z_f = 10 and T_U(0) = 13.67 Gyr and T_G = 13.187 Gyr. Whenever needed we pass from one to the other. In order to minimize the number of free parameters in each model, we assume the following: all galaxies are born at the same redshift z_f = 10; the collapse timescale of τ = 1 Gyr for all galaxies; the Salpeter initial mass function (in number) with a slope x = −2.35; a fraction of stars more massive than 1 M_⊙ equal to ζ = 0.30; absence of galactic winds. However, a few cases are shown for different values of τ, different values of z_f, and in the presence of galactic winds.

4.2. Model results

In this section we discuss the galaxy models obtained with the above prescription for the infall scheme and star formation in particular. First we present the reference case with τ = 1 for all galaxy masses and the prescription for $\stackrel{\sim }{\nu }$ given by Eq. (36), together with the corresponding case $\stackrel{\sim }{\nu }=1$ (which we refer to as the reference case). Then we discuss some cases in which the effect of galactic winds energized by supernova explosions (both Type Ia and Type II) are taken into account. Table 3 lists the models we considered and presents some characteristic features at the last stage with active star formation; this is either the present age for the models without galactic wind or the age at the onset of galactic wind. In the following we mainly present and discuss the models without galactic winds, limiting the discussion of those with galactic winds to some general remarks.

SFR and SF efficiency. In Fig. 16 we show the history of star formation in M_⊙ yr⁻¹ of galaxies with M_B(T_G) equal to 10⁶, 10⁸, 10¹⁰, and 10¹² M_⊙ yr⁻¹ (black, blue, green, and red, respectively) and variable $\stackrel{\sim }{\nu }$ (solid lines) and $\stackrel{\sim }{\nu }=1$ (dashed lines, the reference case). As expected, the SFR starts low, reaches a peak value, and then goes to low values even though it never extinguishes. The peak value is at an age nearly equal to the infall timescale. Models with variable efficiency do not differ from their corresponding reference case with constant $\stackrel{\sim }{\nu }=1$. The reason for this resides in the value of τ. This is shown clearly for the case where τ is allowed to change with galaxy mass. The SFR efficiency $\stackrel{\sim }{\nu }$ is a dimensionless quantity, and therefore is the same for all galaxies; it varies with time from the initial top value 4.04 to 1, as shown in Fig. 17.

Fig. 16. Star formation rate histories of models with MB(TG) equal to 106, 108, 1010, and 1012 M⊙ (black, blue, green, and red, respectively), and with variable $\stackrel{\sim }{\nu }$ (solid lines) and $\stackrel{\sim }{\nu }=1$ (dashed lines, the reference case).

Fig. 17. Temporal variation in SF efficiency of the galaxy models with MB(TG) equal to 106, 108, 1010, and 1012 M⊙ (black, blue, green, and red, respectively). The efficiency is the same for all models, and goes from 4.04 to 1. To better see the curves, the data of each case have been shifted by 0.01 with respect to the others.

The advantage of our choice for the SF efficiency $\stackrel{\sim }{\nu }$ is that this important physical quantity is no longer a free parameter. It is deeply driven by the galaxy mass building process and the timescale associated with it. With our choice for τ, the SF efficiency very quickly reaches its asymptotic value (within about 2 × τ). If τ is increased the timescale over which $\stackrel{\sim }{\nu }$ goes to the asymptotic value gets accordingly longer.

4.2.1. The ratio M_D/M_s and the radius R_e

The stellar radius R_e depends on the dark mass-to-stellar mass ratio M_D/M_s. As already explained this ratio is determined at each time from the current value of the stellar mass built up by star formation, the current mass of baryonic mass M_B(t) and the current mass of dark matter associated with it given by M_D(t) = f_cM_B(t). The ratio M_D/M_s is shown in Fig. 18 as a function of mass M_D (top panel) and age in Gyr (bottom panel). In each galaxy the ratio starts very high and, as time increases, tends to the limit value f_c ≃ 6 as the whole baryonic gas mass is turned into stars by star formation. The general behavior of M_D/M_s as a function of M_D and age is the same to the point that in the bottom panel all the curves overlap. In addition, in this case the ratio M_D/M_s is not an external parameter, but it is determined in a self-consistent way by the internal properties of the models.

Fig. 18. Ratio m = MD/Ms as a function of the mass of dark matter MD (top panel) and age (bottom panel) for our model galaxies. Masses are in M⊙ and ages in Gyr. The model galaxies are in different colors and ranked according to their total baryonic mass MB reached at the present time (black: 106; blue: 108; green: 1010, and red: 1012). In the bottom panel all the lines overlap each other.

With the aid of the m-ratio and the MRR of Eq. (31) we derive the radius R_e of the stellar component M_s and build the MRR of our model galaxies shown in Fig. 19 both along their evolutionary history (black line of filled squares, one for each time step; the present time is at the top and the initial stage at the bottom). Each curve corresponds to a model with a different final total baryonic mass M_B, namely 10⁶ M_⊙, 10⁸ M_⊙, 10¹⁰ M_⊙, 10¹² M_⊙ (from left to right). Although the MRRs of the models are in fair agreement with the bulk of observational data and other theoretical MRRs, a closer inspection of the issue reveals that our theoretical radii are likely overestimated by a factor that is difficult to assess. Our best estimate is about a Δlog R_e ≃ −0.6 to −0.8. The mean radii should be a factor of 4 to 6 smaller. There are many possible causes for this disagreement. First, in addition to the m-ratio in the term (M_D)^1/3, Eq. (31) contains other terms that are affected by some uncertainty. The terms in question are the ratio S_S(n_S)/0.34, the ratio (1.5f_σ)², and finally the ratio (25/m). The first two are simply assumed to be equal to one, while the last contains the ratio m and deserves some remarks. It is clear that it has been introduced as an adjustment factor based on some estimates of the m-ratio derived from current theoretical N-Body Smoothed Particle Hydrodynamic (NBTSPH) simulations of galaxy formation in which only a small fraction of the available gas was used to form stars (see, e.g., Chiosi & Carraro 2002), which explains the factor 25. The present infall models have a different behavior because nearly all the gas is used up to form stars and the limit value of the m-ratio is about f_c ≃ 6. This implies that the above adjustment factor should become (f_c/m), and a consequent reduction in the estimated radius by a factor of about 4 to 6. However, instead of forcing the radius to strictly agree with the data (thus introducing some ad hoc adjustments), we keep the radii as they are, but also keep in mind that in reality they could be 4 to 6 times smaller than estimated. This would immediately affect our evaluation of the specific intensity $I_{\text{e}}=L/(2\pi R_{\text{e}}^2)$ that could be a factor of 16 to 36 higher than our straight evaluation (see below).

Fig. 19. Mass-radius relations (MRRs) of our model galaxies labeled by their present-day total baryonic mass MB(TG) equal to 106, 108, 1010, and 1012 M⊙ (from left to right). Each line made of filled black squares represents the whole evolutionary history of Ms and Re, both increasing with time (present-day at the top). These models are compared both with observational and theoretical data from different sources: (i) the observational data of Burstein et al. (1997) from GCs to GCGs (light blue dots) and the ETGs by Bernardi et al. (2010, red dots); (ii) the Illustris-1 galaxies (light green dots); (iii) the low initial density models (blue squares and their best fit) and the high initial density models (red squares and their best fit) by Chiosi & Carraro (2002; iv) the early hierarchical models by Merlin et al. (2012, black squares and their best fit); (v) the Fan et al. (2010) MRRs for different values of the formation redshift (from top to bottom) zf = 0, 1, 5, 10, and 20; and (vi) the MRR by Chiosi et al. (2020, dark yellow line). See the text for more details.

4.2.2. Other important relationships: L versus R_e and R_e versus σ

The uncertainty on the radius also affects other important relationships, such as the luminosity-radius relation (LRR) shown in Fig. 20 and the radius velocity dispersion relation (RSR) displayed in Fig. 21. The theoretical and observational data are compared by Burstein et al. (1997) and Bernardi et al. (2010); for these latter the mean color (B − V) = 0.85 was applied to the M_V magnitudes to get the B-luminosity. In both cases the best results are for radii reduced by a factor of 4.

Fig. 20. B-luminosity–radius relation (LBRRs) of our model galaxies with MB(TG) equal to 106, 108, 1010, and 1012 M⊙ (from bottom to top). For each mass two lines are shown: one with the original radii (dashed black line) and the case with the radii decreased by a factor of 4 (line made of filled black squares; see text). The models are compared both with observational data from Burstein et al. (1997) from GCs (magenta dots) to Dwarf Galaxies (blue dots), to ETGs (red dots), GCGs (light blue dots), and the ETGs by Bernardi et al. (2010, red dots, which overlap the previous ones). The agreement for the smaller radii case is soon evident (see the text for more details).

Fig. 21. Re–σ relation (RSR) of our model galaxies. In this figure the same models, observational data, color-coding, and symbols as used in Fig. 20 are adopted.

4.2.3. Remarks on the luminosity

Before proceeding further it is worth commenting on the luminosity of the model galaxies. As already explained, for the sake of a quick assessment of the model galaxies luminosity in the B and V passbands, we used suitable linear relationships between the absolute B and/or V magnitudes of the Johnson system and the mean age T based on SSPs of mean metallicity. One may argue that the luminosities derived in this way are much different from those evaluated by means of full population synthesis technique (i.e., by integrating the spectral energy distribution of SSPs over the SFR, initial mass, and the metallicity range spanned by the stellar populations of galaxies at each time; see Bressan et al. 1994, for all details). This is done a posteriori, once the whole SFR(t), M_s(t), and metallicity Z(t) are known. The results are shown in Fig. 22 for the M_B(T_G) = 10⁶, 10¹² M_⊙, where the solid lines are the luminosity from the analytical relationships, and the dotted lines the luminosity from full population synthesis. The two luminosities differ from each other by a maximum factor of 3 in the far past when the SFR(t) was highest (ages of about 1.5 Gyr), while they coincide in the less remote past (roughly the past 5−6 Gyr). Therefore, our approximation that nicely speeds up the model calculation is reasonable and leads to acceptable results. Our luminosities can be safely used for the present purposes.

Fig. 22. Luminosity LV vs. age (in Gyr) relation (LAR) of our model galaxies. The cases MB(TG) = 106 (bottom) and 1012 M⊙ (top) are shown. The solid lines are luminosities derived from the analytical relationships, while the dotted lines are those from full population synthesis.

4.2.4. The I_e versus R_e plane

Together with the FP and the luminosity-velocity dispersion relation, otherwise known as the Faber–Jackson (FJ) relation, the I_e − R_e plane is one of the most studied projections of the FP. The uncertainty on the radius R_e (a factor of 4) reflects on I_e as an increase of a factor of 16 at fixed stellar mass of the galaxy. The results for our models (with no galactic winds) are shown in Fig. 23 and are compared with the observational data of Burstein et al. (1997) from GCs to GCGs using the same color-coding as in Fig. 20. The evolutionary sequences on display are for model galaxies with M_B(T_G) equal to 10⁶, 10⁸, 10¹⁰, and 10¹² M_⊙ (from left to right). For each mass we display two lines, the one with the original radii (dashed black line) and the case with the radii decreased by a factor of 4 and the specific intensity I_e increased by a factor of 16 (line made of filled black squares). The time evolution goes from the top to the bottom of each line. The present-day stage is the bottom point of each line. Finally, the thick dashed line is the border of the ZoE. We note that no model at the present time falls in the ZoE, but all are below it.

Fig. 23. Ie − Re plane of our model galaxies compared with the observational data of Burstein et al. (1997). The color-coding of the data is the same as in previous figures. There are two groups of models. The black thin dashed lines are models with original radii, while the thick lines made of filled black squares are those with the radius decreased by a factor of 4 and the specific intensity Ie increased by a factor of 16. The galaxy mass is MB(TG) equal to 106, 108, 1010, and 1012 M⊙ (from left to right). Along each line the time runs from zero to present age (from top to bottom). The formation redshift of all the models is zfor = 10.

The present models cannot account for the data of GCs (as expected). Even if the model with M_B(T_G) = 10⁶ M_⊙ crosses the region of GCs, it cannot reproduce these objects because the present radius and specific intensity are too large and too low, respectively. The ongoing star formation yields objects that are too luminous and too large, and they do not match with general properties of GCs. What is needed are models in which star formation ceased and the radius stopped growing long ago (short initial episode followed by quiescence perhaps because of strong galactic wind), or to take into account the important transformations induced by the interaction with the Galaxy.

Similar considerations can be applied to clusters of galaxies for which different types of models should be set up. To develop a suitable model for the formation and evolution of GCGs along the same lines we followed for single galaxies is beyond the aims of this study and we leave the subject to a future investigation.

The β − $L_0^{\prime }$ space. With the aid of Eqs. (3)–(14) we derive the exponent β and proportionality factor $L_0^{\prime }$ along the whole evolutionary sequence of our model galaxies evolved without galactic winds. In the left panel of Fig. 24 we present all cases under consideration: models with larger radii and models with smaller radii (factor of 4) for the two photometric passbands in use (B and V Johnson). Along each line the time increases from the bottom to the top, and the last stage at the present age is indicated by the mass label (total asymptotic baryon mass). For each galaxy mass (M_B(t_G)) the results are nearly the same; all sequences overlap each other (see also the entries of Table 4 containing the slope α and zero-point γ of their linear best fits). It turns out that the relationships in question depend only on the galaxy mass. Remarkably, the exponent of the $L=L_0^{\prime }{\sigma }^{\beta }$ relation (β) is positive during the early stages and negative afterwards. The luminosity first increases with σ and then decreases with it afterwards. Finally, we note that all curves cross each other at β ≃ 3 and log $L_0^{\prime }$ ≃ 2.5, values very close to the observed FJ relation. To confirm this, in the right panel of Fig. 24 we plot the same relations for the Illustris-1 models grouped at different redshifts from z = 4 to z = 0. Now the situation is not the same as before; in each group with the same redshift, mass and age vary from galaxy to galaxy. Furthermore not all masses are present at each redshift: samples at high redshifts (e.g., ≥1,  6) are dominated by low mass objects (masses lower than 10⁸ M_⊙ are missing anyway because of the mass resolution); massive objects up to 10¹² M_⊙ are present at lower redshifts. However, the resulting distributions in the β − $L_0^{\prime }$ plane are very similar to that of the left panel, and remarkably there is also some evidence of the β ≃ 3 cross-point. This strongly supports the notion that infall models nicely mimic the numerical hierarchical simulations.

Fig. 24. The mutual relationship between β and $L_0^{\prime }$. Left panel: $L_0^{\prime }$ − β relation of our models (left panel) evolved without galactic winds. All the relationships are nearly identical passing from models with large radii to those with smaller radii (by a factor of 4), and the corresponding solutions of the equations for β and $L_0^{\prime }$, and finally changing only the photometric passband in use. The relationships seem to depend only on the galaxy mass. Right panel: $L_0^{\prime }$ − β relation for the artificial galaxies of Illustris-1. Each color corresponds to a different redshift epoch: green (z = 4), blue (z = 3), yellow (z = 2.2), brown (z = 1.6), magenta (z = 1.0), dark gray (z = 0.6), red (z = 0.2), and light gray (z = 0).

The most important relation to look at and to examine in detail is luminosity versus velocity dispersion. This is shown in Fig. 25 which displays log(L_B/L_⊙) versus log σ for the model galaxies and compares it with observational data of Burstein et al. (1997). In the main panel we show three possible relationships: (1) The first is the plain L_B/L_⊙ versus σ of the models with their original luminosity and radii (thick black curves). Along each curve the evolution starts at the bottom point of each line and proceeds to the final stage indicated by the label M_B(t_G) in solar units. We note that during the galaxy lifetime the L_B/L_⊙ versus σ relation bends over past a certain age toward lower luminosities and lower velocity dispersion. This happens roughly past the peak of star formation. While the luminosity decrease can be easily understood, the decrease in velocity dispersion of the stellar component needs some explanation. Stars during their lifetime can explode as Type II and Type I supernovae. In the first case a small remnant is left (neutron star or black hole), in the second no remnant at all is left. They can also lose a great deal of mass by stellar winds. In any case the total mass in stars is expected to decrease, and so does the velocity dispersion; (2) The second case is the associated $L_B/L_{\odot }=L_0^{\prime }{\sigma }^{\beta }$ relation in which the original β and $L_0^{\prime }$ are used (red curves) together with the linear fit limited to the descending branch of each curve (red solid lines); (3) Finally, the $L_B/L_{\odot }=L_0^{\prime }{\sigma }^{\beta }$ relation, in which the correction on the radius is applied and new values of β and $L_0^{\prime }$ are derived. It is worth recalling that the $L_0^{\prime }$ versus β relation remains unchanged. The results are shown by the green curves. The small inset in Fig. 25 shows the case of the M_B(t_G) = 10¹⁰ M_⊙ in more detail for the sake of clarity. Similar results are found for the V passband that are not shown here. Two important features are soon evident. First, the relations log(L_B/L_⊙) versus log σ, based on the model history past the star formation activity period, have a similar slope but different zero point (which depends on the galaxy mass). The manifold of these relations provides a sort of natural width to the luminosity–sigma relationship. The mean slope of the manifold agrees with the current value of the observed FJ. Second, the theoretical relations marginally agree with the body of ETGs; a steeper slope at luminosities above log L_B/L_⊙ ≃ 9 would be more appropriate.

Fig. 25. Luminosity LB/L⊙ vs. velocity dispersion σ (in km s−1) relation of our model galaxies. For each mass (labeled MB(tG), as indicated) three relations are shown: (1) original models with no revision of the radii (lines made of filled black squares); (2) models whose luminosity is derived from the $L_B/L_{\odot }=L_0^{\prime }{\sigma }^{\beta }$ relation with the original β and $L_0^{\prime }$ (curves made of red squares) together with the linear fit limited to the descending branch of each curve (black solid lines); (3) models in which the radii have been revised and new values of β and $L_0^{\prime }$ are calculated (green curves). In each case the luminosity–sigma relation bends past the stage that roughly corresponds to the maximum stellar activity. From this stage the luminosity and velocity dispersion decrease (see text for details). The inset shows the case of MB(tG) = 1010 M⊙ for the sake of clarity. The models are compared with the data by Burstein et al. (1997) from GCs to GCGs; the same color-coding as in the previous figures is used.

The simplicity of the current models cannot lead to better results. A possible improvement could be given by allowing small secondary episodes of star formation. The argument is as follows. The luminosity is the product of the star mass times the flux per unit mass. We call L_o = f_oM_o the original luminosity and L_n = f_nM_n the expected luminosity including some recent star formation activity; L_nM_n is in turn f_yM_y + f_oM_o, where f_yM_y is the contribution by the episodic stellar activity; it follows that (f_yM_y + f_oM_o)/(f_oM_o) = λ = L_n/L_o. Based on the current observations we would expect approximately λ ≃ 2. Now we may also assume M_y ≪ M_o so that the total stellar mass and velocity dispersion in turn remain nearly constant. Indicating with θ = M_y/M_o, we gets f_y/f_o = λ − 1/θ ≃ 5 − 10, which is not impossible according to current population theories for SSPs. Third, in the theoretical models the exponent β of the $L_B/L_{\odot }=L_0^{\prime }{\sigma }^{\beta }$ relationship (a generalization of the FJ) can be either positive or negative depending on the particular evolutionary stage of the galaxy. Therefore, among the observational data both values of β are to be expected without violating the trend indicated by the FJ (i.e., that the luminosity of galaxies increases with the velocity dispersion, hence the mass of the galaxy).

Galactic winds. Long ago, Larson (1974) postulated that the present-day color–magnitude relation (CMR) of ETGs could be the result of galactic winds powered by supernova explosions, and thus initiated a long series of chemo-spectro-photometric models of elliptical galaxies based on this idea (see, e.g., Tantalo et al. 1998a, and references). In brief, gas is let to escape from the galaxy, and star formation is supposed to halt when the total thermal energy of the gas equates its gravitational binding energy. This idea has been extended including the effect of stellar winds in the thermal energy budget of the gas. It was also included in NBTSPH models of galaxies (see Merlin et al. 2012, and references).

The same scheme proposed by Tantalo et al. (1998a) is adopted here, but with minor modifications because of the much simpler present formalism. As already said, the thermal energy of the gas is the sum of three contributions, namely type I and II supernovae and stellar winds from massive stars,

$$\begin{array}{c}\hfill E_{\mathrm{th}}(t)=E_{\mathrm{th}}{(t)}_{\mathrm{SNI}}+E_{\mathrm{th}}{(t)}_{\mathrm{SNII}}+E_{\mathrm{th}}{(t)}_W,\end{array}$$(39)

where each term has the generic expression

$$\begin{array}{c}\hfill E_{\mathrm{th}}{(t)}_j={\displaystyle {\int }_0^t}{ϵ}_j(t-t\prime )R_j(t\prime )M_{\mathrm{B}}(t_{\mathrm{G}})\mathrm{d}t\prime \end{array}$$(40)

with j = SNI, SNII, W (using the obvious meanings for the symbols). The normalization factor M_B(t_G) in these equations is required to calculate the energy in physical units. The time t′ is either the SN explosion time or the time of ejection of the stellar winds, as appropriate. The functions ϵ_SN(t) and ϵ_W(t) are cooling laws governing the energy content of supernova remnants and stellar winds, respectively. Finally, star formation and chemical enrichment are halted, and the remaining gas content is supposed to be expelled out of the galaxy (winds) when the condition

$$\begin{array}{c}\hfill E_{\mathrm{th}}(t)\ge {\mathrm{\Omega }}_{\mathrm{g}}(t)\end{array}$$(41)

is verified. For all other details concerning the above rates, the evolution of SN remnants and stellar winds, and how much of the initial energy budget is shared with the gas to energize it, and finally for the expression for the gravitational energy of the gas in presence of baryonic and dark mass and their space distribution in a galaxy, we refer to Tantalo et al. (1998a).

A small sample of models with galactic winds are calculated and their main features are summarized in Table 3. It is worth noting that the onset of galactic winds occurs at younger and younger ages as the galaxy mass increases. Thanks to this effect, the models obey the constraint imposed by observational data on chemical elements like carbon (C), oxygen (O), and magnesium (Mg), also known as α-elements, and iron (Fe) and their ratios [α/Fe]. The high mass galaxies are more α-enhanced ([α/Fe] > 0) than the low mass ones ([α/Fe] ≤ 0). This cannot be easily reconciled with other properties of the same objects (see Chiosi et al. 1998; Tantalo & Chiosi 2002 for detailed discussions of this issue and possible ways out). In the present models we take the suggestions by Chiosi et al. (1998) and Tantalo & Chiosi (2002) into account.

4.2.5. The role of galactic winds

The main lines of the discussion for models without galactic winds holds also for the new models. Therefore, we focus on key relations such as the I_e − R_e plane, which is shown in Fig. 26. A comparison with the same plot in Fig. 23 shows that there is no visible difference when passing from models without galactic winds to those with. The reason for this is the kind of star formation at work. Because of the short infall timescale and the dependence of $\stackrel{\sim }{\nu }$ on the inverse of the velocity dispersion, most of the stars are in place before the occurrence of galactic winds. To somewhat alter this trend the parameter τ should be changed and made to depend on the galaxy mass, for instance long in low mass galaxies and short in the high mass ones. To further investigate this point is beyond the aims of this study.

Fig. 26. Ie − Re plane of our model galaxies with galactic winds powered by the energy input from supernova explosions and stellar winds. There is no visible difference with respect to the same plane of models without galactic wind. The same notation, symbols, and color-coding as in Fig. 23 is adopted.

Role of the initial mass function. To avoid misunderstanding, we recall here that the present models are calculated with the classical IMF of Salpeter (1955). Therefore, the M_s/L ratio based on these models has this fundamental limitation and cannot by applied to investigate the problem of the FP tilt in a very general way. Our infall models can easily be adapted to include the popular IMFs in the literature different from the Salpeter case (see, e.g., Chiosi et al. 1998), where the IMF is let vary with the physical condition (mean density, temperature, and velocity dispersion) of the gas inside a galaxy, and therefore with time for a galaxy of given mass, and with time and mass in objects of different mass. However, for the aims of this study, in order to simplify the model, we thought it wise to rely on the classical IMF of Salpeter. If the present models were applied to the issue of the FP tilt, most likely they could account for only half of the observational tilt. This subject was specifically addressed in Chiosi et al. (1998) with good results for the tilt of the FP of ETGs in the Virgo and Coma clusters.

4.2.6. Changing the galaxy mass and z_f

An important feature of the models is related to the formalism in use. According to the formalism and equations widely described in Tantalo et al. (1998a) all relevant physical quantities describing the model and its temporal evolution are suitably normalized to the asymptotic baryonic mass M_B(t_G), for instance the gas mass at time t is expressed as G_g(t) = M_g(t)/M_B(t_G), equally for the star mass G_s(t) = M_s(t)/M_B(t_G), and the current total baryonic mass G_B(t) = M_B(t)/M_B(t_G). The amount of dark matter at any time is simply related to the current baryonic mass via the cosmic ratio; the components are mixed together so that they fall together at the same rate. Furthermore, the accretion rate and the SFR, for example, are all expressed using the same kind of normalization. The advantage is that the timescale of mass accretion τ, the cosmic ratio f_c, and the rest is parameter-free, so that the only free quantity is the asymptotic baryonic mass M_B(t_G). This allows us to generate models for any value of M_B(t_G).

All galaxy models discussed so far are calculated assuming the redshift of galaxy formation z_f = 10, although other values are possible. Higher values are unlikely, whereas lower values are more plausible. Since the age of the Universe T_U depends only on the cosmological model in use and is therefore a fixed quantity, the age of a galaxy T_G expressed by T_G = T_U − T_U(z_f) at decreasing z_f becomes younger. Consequently some features of the models will change, such as total ages, radii, luminosities, and specific intensities. The following values of z_f are considered: 5, 3, 1, and 0.5, in addition to the previous set with z_f = 10. The results are shown in Fig. 27 limited to the cases z_f = 5 (black lines) and z_f = 0.5 (red lines). The case z_f = 10 runs nearly over the case z_f = 5. All the others are between the case z_f = 5 and z_f = 0.5. From left to right, the galaxy mass is M_B(T_G) = 10⁶10⁸, 10¹⁰, 10¹² M_⊙. The age increases along each line from the top to the bottom. The final age (in Gyr) decreases from 13.19 for z_f = 10 to 12.47 for z_f = 5, 11.48 for z_f = 3, 7.73 for z_f = 1, and 5.02 for z_f = 0.5 (see Table 5 for more information on the final stage of each model). In Fig. 27 the final stages are represented by the green circles (some of which overlap).

Fig. 27. The Ie − Re plane of our model galaxies with different formation redshift zf, namely 10, 5, 3, 1 and 0.5. The four green points of different colors are the present-day stage of model galaxies whose existence began at redshifts from 0.5 to 5. The effect is quite small.

From these data we derive that variations in z_f from 10 to 0.5 yields variations in log(I_e) of about Δlog(I_e)≃0.5, while the radius does not change significantly. More efficient star formation in recent times generates more luminosity, and hence higher specific intensity I_e. This is achieved by changing τ from 1 to 5 Gyr (in the case of the 10¹⁰ M_⊙ galaxy) yielding Δlog(I_e)≃0.4. Recent bursts of star formation either by internal causes or mergers would also increase I_e. Analyzing all the implications of this is beyond the aims of this study. What we can say with confidence is that a significant scatter in the I_e − R_e plane is likely to occur. In any case, the gross distribution of galaxies in this plane (but for GCs and GCGs) is accounted for by these models. Finally, the homologous behavior of the models and the limited effect of the formation redshift on their evolutionary behavior make it possible to generate simulations of the distribution of the large number of galaxies in the parameter space we are investigating, in practice at no cost.

4.2.7. A test of consistency

The galaxy models we have presented are based on physical assumptions such as the infall picture, the SFR, the mass–radius relationship, and the population synthesis governing their luminosity in different passbands, which are not explicitly related to our interpretation of the parameter space of galaxies (luminosity, stellar mass and radius, velocity dispersion, and specific intensity), the FP in multi-dimensional space, and its possible projections onto different planes that led us to the $L=L_0^{\prime }{\sigma }^{\beta }$ relationship with $L_0^{\prime }$ and β changing from galaxy to galaxy and for each of them also with time. In light of this, we have made some detailed predictions about $L_0^{\prime }$ and β and derived a number of equations whose solutions on the one hand yield $L_0^{\prime }$ and β as functions of L, M_s, R_e, σ, and on the other hand allows us to construct the expected relationships among pairs of fundamental variables, such as I_e − R_e, I_e − σ, and R_e − σ. Among these we choose here as an example the variables I_e and R_e and compare the values given by the models with those derived from Eqs. (16) and (17). The comparison is shown in Fig. 28; the galaxy mass M_B(T_G) is 10⁶, 10⁷, 10⁸, 10¹⁰, and 10¹² M_⊙ (the case M_B(T_G) = 10⁷ M_⊙ is added). On the abscissa are the input values from the models (labeled I_e[i] and R_e[i]) and on the ordinate the values calculated from Eqs. (16) and (17), (labeled I_e[c] and R_e[c]). In general there is a surprisingly good agreement between the [i] and [c] quantities, but for some particular stages at which the [c]-values rapidly diverge and change sign. The cause lies in the analytical relationships themselves that contain various exponents (e.g., γ, [(β − 2)/(1 + 3/γ)], [(β − 2)/(3 + γ)]), which in turn are functions of β, which varies in the course of evolution. In this narrow interval the disagreement is of mathematical nature with no physical implications. It simply means that these analytical relationships cannot be safely used to derive the corresponding variables.

Fig. 28. Comparison of Ie and Re derived for the model galaxies (indicated by the suffix [i] and those calculated with relations (16) and (17) for galaxies with asymptotic baryonic mass MB(TG) = 106, 107, 108, 1010, and 1012 M⊙ from left to right. The redshift of galaxy formation is zf = 10.

4.2.8. General remarks and preliminary conclusions

Since the [i]- and [c]-values are nearly coincident, using the analytical relationships would predict results in the various projection planes we have examined that are identical to those obtained from using the numerical galaxy model. The overall agreement between the model and analytical approach lends strong support to the idea at the base of the analytical view, where the relation between the luminosity and velocity dispersion of a galaxy is governed by $L=L_0^{\prime }{\sigma }^{\beta }$ in which both β and $L_0^{\prime }$ vary with the galaxy mass, evolutionary stage, and hence time and redshift, and that these quantities in turn are intimately related to key physical parameters such as the stellar mass and radius, the velocity dispersion (a measure of the gravitational potential well), the SFR, the infall timescale, and finally the ratio M_D/M_s. The distribution of galaxies on the usual diagnostic planes (e.g., FP, FJ, I_e − R_e, M_s − R_e, L − R_e, R_e − σ, and the border of the exclusion zone) mirrors the mean behavior of galaxies, each of which has its particular history and is observed at some evolutionary stage.

5. The important role of β

In this section we cast light on the role of β. To this end we adopt the reference case (z_f = 10 and τ = 1) and leave the issue of galactic winds aside. For this case we present a few basic relationships among β and other important parameters, namely the SFR (in M_⊙ yr⁻¹), the age (in Gyr), and the specific intensity I_eB or I_eV (in L_⊙ pc⁻²). These relationships are shown in Fig. 29. In the left panel the homologous nature of the galaxy models is evident. All curves have the same shape, but each one is separated from the others by the homology parameter, namely the total baryonic mass at the present age M_B(T_G) annotated along each curve. The temporal evolution occurs from the top to the bottom of each curve. Identical behaviors are found between β and the luminosity L_B or L_V (in L_⊙), and the velocity dispersion σ (in km s⁻¹). However, these relationships are not shown here for the sake of brevity. The central panel of Fig. 29, showing the variation in β with age, still displays the dependence of the results on the homology parameter, and thus there are four different curves, one for each value of M_B(T_G). Finally, in the right panel we show the dependence of β on the surface brightness I_eB; all curves collapse to a single relation, and the homology is destroyed by the underlying relationship between the mass and the effective radius of the models. A similar relation is found between between β and I_eV. The analytical relations between β and I_e are given by

$$\begin{array}{cc}\hfill \beta & =3.159\mathrm{Log}(I_{\mathit{eB}})-2.003,\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill \beta & =3.159\mathrm{Log}(I_{\mathit{eV}})-1.900.\hfill \end{array}$$(43)

Fig. 29. Correlations of β with SFR, Age and Ie. Left panel: relationship between β and SFR. Each curve labeled MB(TG) is identified according to the color-coding adopted in the previous figures. The total baryonic mass is the homology parameter separating each curve from the others. More similar trends are found for the luminosity LB and LV, and for the velocity dispersion σ, which are not shown here for the sake of brevity. In all three relations β mirrors the behavior of the SFR, the luminosity in turn, and finally the velocity dispersion. The SFR is in M⊙ yr−1. Middle panel: relationship between β and age (in Gyr). Symbols and color-coding have the same meaning as in the left panel. Right panel: relationship between β and IeB in L⊙ pc−2. The lines corresponding to different masses of galaxies have been shifted by 0.1; they collapse to a unique curve given by β = 3.159Log(IeB)−2.003. The long dashed line is the best fit of the theoretical data. An identical relation can be found for IeV with the same slope but slightly different zero point.

The evolution along each line is from top right to bottom left, and the present stage is the last point where M_B(T_G) is annotated. The above relations indicate both the path followed by a single galaxy in the course of its history and also the locus on the β − I_e plane of galaxies of different mass observed at the present age. There is no appreciable effect of different formation redshifts, at least in the interval 10 ≥ z_f ≥ 1, nor of different accretion timescale τ in the interval 1 ≤ τ ≤ 5 Gyr. We estimate a total effect on β by redshift z_f and accretion timescale τ on the order of Δβ ≃ 2 over the interval of I_eV of interest here.

The linear relation between β and I_e shown by our models is a very intriguing result that demands a thorough analysis because observational data and numerical hierarchical models seem to indicate a different picture. The situation is best illustrated by Fig. 30 comparing data and models from different sources. On the observational side we have three data sets: Burstein et al. (1997), WINGS, and Bernardi et al. (2010). The last two (mainly devoted to ETGs) are based on equivalent methods to estimate R_e, and therefore yield similar results for the β − I_e plane. In contrast, the first one contains objects going from GCs to DGs, LTGs, ETGs, and GCGs, and differs in the method used to derive the effective radius, and consequently yields different relationships in the β − I_e plane. Owing to this, some preliminary remarks are needed. First of all, the data of Burstein et al. (1997) are in the B-band so that must be transformed into the V-band. This is made by means of the relation

$$\begin{array}{c}\hfill log{L}_V=0.4[{(B-V)}_0-0.65]+log{L}_B,\end{array}$$

Fig. 30. Data and theoretical models for the β − IeV plane. Data from the different sources are plotted: (i) Burstein et al. (1997, light blue). Three sequences are seen: GCGs and GCs, ETGs (no evidence of star formation), LTGs and DGs (evidence of ongoing star formation). By construction, these data are well behaved with no evidence of dispersion. (ii) WINGS data (red squares) showing large dispersion in both coordinates, log Ie is always positive, and β can be very large both positive and negative. (iii) ETGs of Bernardi et al. (2010, open green circles), but limited to z ≤ 0.02. The Illustris-1 model galaxies are indicated by the blue dots; their distribution closely mimics that of the observational data. Finally, the long dashed red line shows the present-day position on the β − Ie plane of our models for the reference case (with τ = 1 Gyr, zf = 10, and no galactic winds). This line coincides with the lower border of the Illustris-1 distribution in the β > 0 hemi-plane.

where the luminosities are in solar units, (B − V)₀ is the color, and −0.65 is the difference between the B and V photometric constants (5.48 and 4.83, respectively). Second, recalling that the luminosity L_e given by Burstein et al. (1997) is the amount of light falling within the effective radius R_e, where half the total luminosity is found, we scale it by a factor of two to make it consistent with the definition of I_e we have adopted.

The observational data for Re, L_V, M_s, (M_s/L_V), σ, and I_eV are fed into the system of Eq. (13) and the solutions β and $L_0^{\prime }$ are derived. The light blue points in Fig. 30 are the Burstein et al. (1997) data; three sequences are seen: the GCG-GC sequence, the one of ETGs (no evidence of star formation), and the one of LTGs and DGs (evidence of ongoing star formation). By construction, the data of Burstein et al. (1997) are well behaved with no evidence of dispersion. The red squares are the WINGS data showing large dispersion in both coordinates, log I_e is always positive, and β can be very large both positive and negative. Very similar results are found with the Bernardi et al. (2010) data, the open green circles. The Illustris-1 model galaxies are indicated by the blue dots; their distribution closely mimics that of the observational data. Finally, the long dashed red line shows the present-day position on the β − I_e plane of our models for the reference case (with τ = 1 Gyr, z_f = 10 and no galactic winds). This line coincides with the lower border of the Illustris-1 distribution in the β > 0 hemi-plane. Choosing different values of z_f in the interval 0.5 ≤ z_f ≤ 10 does not significantly shift the line predicted at the present time; the same is true for the effect of galactic winds. Lumping all these effects together we expect a typical width of this border line of about Δβ ≃ 10 over the I_eV interval of interest here.

From this preliminary comparison we can conclude that there is mutual consistency among the different sources of data and of these latter with models. The major issue now is to understand the physical causes of the large dispersion in β for all values of I_e ≥ 1. Looking at Eq. (14), providing the solutions β and log($L_0^{\prime }$) of our equations, we note that under suitable conditions the term 1 − 2A′/A at the denominator of log($L_0^{\prime }$) gets very close to zero, consequently β can be either very large and positive or large and negative. As already discussed in Sect. 3, when this happens the system is under conditions of strict virialization. The sign of β depends on the particular history of the constituent variables (M_s, R_e, L, and I_e), in other words whether the term 2A′/A is tending to 1 from below (β > 0) or above 1 (β < 0). From an operative point of view we can define the “state close to strict virialization” when |β|> 20. This would account for the gap in the negative hemi-plane of Fig. 30.

The question arises of whether data and models ever reach the condition of full virialization indicated β ⇒ ±∞ or if they remain somewhat far from it. The answer is that both possibilities occur. On the observational side, given any galaxy for which the set of parameters (M_s, R_e, L, I_e, and σ) has been measured, it is not granted that they would satisfy the virialization condition. The major uncertainties are with M_s and R_e, and in turn I_e. Therefore, many of them crowd in the interval −1 ≤ β ≤ 20, which implies deviations from virial equilibrium. However, with the present data it is not possible to say whether this is due to insufficient accuracy in the parameter determinations or to real deviations from virial equilibrium, for example caused by recent mergers, harassment, loss of mass, interactions. However, there are many other galaxies in both hemi-planes with |β|> 20, which is a strong indication that they are close to virial equilibrium.

On the theoretical side, our model galaxies with infall (without dynamics) seem to be in a state far from strict virialization. This is suggested by the small values of β reached at the present time. The reason for that lies in the way the models are built up. In brief, a mass point description with no dynamics is adopted, the total mass is assigned (via the accretion law), the stellar mass is derived from star formation, the effective radius is estimated from a suitable relationship, the luminosity is evaluated from the stellar mass, and a mean luminosity–age relationship for a fictitious SSP with mean metal content is derived (the difference with respect to the luminosity correctly derived from the theory of population synthesis via the history of star formation and metal enrichment is not large but still significant). Finally, the velocity dispersion is derived from the VT with the current values of M_s and R_e. The major uncertainties are in R_e and L. Therefore, our set of basic parameters is not necessarily able to fulfill all the requirements imposed by the VT. In consequence, our β values are always quite small (e.g., smaller than 25−30), implying that full virialization is not reached. However, this failure is not as severe as it appears because small adjustments of R_e and L are possible, while the models are successful in many other aspects.

The situation is much better with the Illustris-1 models. If a good number of galaxies have β < 25 − 30, as is the case of our models, a large number of objects is still clearly seen in the regime of strict virialization because of their high positive and/or negative βs. The inclusion of real dynamics and the hierarchical scenario at work provide much better conditions to bring the action of virialization into evidence. In the hierarchical scenario, mergers, ablation of stars and gas, harassment, secondary star formation, and inflation of dimension by energy injections of various kinds induce strong variations in the structural parameters and hence strong temporary deviations from the virial conditions. However, once this has happened the viral conditions can be soon recovered over a suitable timescale. This can be short or long depending on the amount of mass engaged in the secondary star formation activity and the amount of time elapsed since the star formation event took place (see the burst experiments in Chiosi & Carraro 2002; Tantalo & Chiosi 2004b). As a consequence, detecting systems on their way back to virial equilibrium is likely a frequent event, and thus explains the high dispersion seen on the β − I_e plane.

In principle, the value of β evaluated for each galaxy could provide a useful hint about the equilibrium state reached by the system. Most likely, the condition of strict virial equilibrium is a transient phenomenon that could occur several times during the life of a galaxy. This is perhaps suggested by the large numbers of galaxies with both low and positive values of β and high positive and/or negative values of β.

6. Discussion and conclusions

The aim of this paper is to prove that the difficulties encountered in understanding the distribution of galaxies on the FP in the parameter space σ, R_e I_e, and its projections on the three coordinate planes, can be removed by introducing the $L=L_0^{\prime }{\sigma }^{\beta }$ relation as a proxy of evolution, in which β and $L_0^{\prime }$ vary from galaxy to galaxy and for each of them in the course of time (see D’Onofrio et al. 2017; D’Onofrio et al. 2019; D’Onofrio et al. 2020; D’Onofrio & Chiosi 2021, for previous efforts along this line of thought). The continuous variation in β and $L_0^{\prime }$ traces the path followed by each ETG in the L − σ plane. The $L=L_0^{\prime }{\sigma }^{\beta }$ law together with the VT yield a set of relations R_e − σ, I_e − σ, and R_e − M_s that nicely reproduce the data and suggest the existence of a system of two equations in the unknowns β and $L_0^{\prime }$ with coefficient functions of M_s, R_e, L, and I_e that for each galaxy determine the value of β and $L_0^{\prime }$. With the aid of these relations we can determine the instantaneous position and direction of a galaxy on the FP and the projection planes.

The analysis is made in two steps. In the first step the problem is addressed from an observational point of view, inferring from the data the expected position and evolutionary direction of a galaxy in the various planes and owing to the large number of galaxies in the samples, the range of values spanned by β and $L_0^{\prime }$ is determined. In the second step simple models of galaxy formation, structure, and evolution are set up, and the basic equations are solved at each time step of a galaxy’s lifetime so that the history of β and $L_0^{\prime }$ is known. Based on these results, the various projection planes are examined and consistency between observational data and theoretical models is found. The same procedure is applied to the literature galaxy models calculated in the framework of the hierarchical scenario. The theoretical results are compared with the observational data and good mutual agreement is found. Based on this, we conclude that the starting hypothesis about the real existence of the $L=L_0^{\prime }{\sigma }^{\beta }$ relation is correct. In more detail, the present analysis has clarified the following issues:

1. The FP can be understood as the average of the single FP-like relations valid for each galaxy (Eq. (4)). The coefficients of the FP-like relation are functions of β. This means that the FP must evolve with redshift and that its coefficients depend on the adopted waveband (in which observations are taken) and on the nature of the data sample (how many ETGs are included) as confirmed by the current observational data.

2. All the features of the FP projections can be explained at the same time. This includes the curvature of the relations, which turns out to depend on the existence of positive and negative values of β, and the existence of the ZoE (i.e., the line marking the separation between the permitted and forbidden regions in these planes). No galaxies can reside in the ZoE.

3. The FP and all its projections, such as the classical Faber–Jackson relation, are instantaneous pictures of the present-day situation. They should change with redshift, and hence the lifetime of the galaxies.

4. The ZoE is obtained in a natural way as the only possible evolutionary path for objects with large positive and negative β values that are well virialized. These objects in general stopped their star formation long ago, and their luminosity is progressively decreasing. When ETGs become passive quenched objects, with a luminosity decreasing at nearly constant σ, the galaxies can only move in one direction, that given by the large positive and negative values of β.

5. The infall model galaxies built here, although lacking the dynamical component, are in reasonable agreement with observations and support the idea that β and $L_0^{\prime }$ vary across time, and therefore that the $L=L_0^{\prime }{\sigma }^{\beta }$ law is a plausible empirical relation accounting for the variation that occurs in a galaxy.

6. All the diagrams built using the structural parameters are sensitive to the temporal evolution of galaxies, simply because each individual object moves in a different way according to the value of β.

7. Observations and theory both suggest that the $L=L_0^{\prime }{\sigma }^{\beta }$ relation provides an empirical way of capturing the temporal evolution of ETGs (and probably of late-type objects) because the values of β are related to the history of mass assembly and luminosity evolution. We are thus tempted to suggest that Eq. (13) are two important equations governing the evolution of ETGs.

8. Finally, the large negative and positive values of β of some galaxies can be considered as the signature that these system are very close to the virial equilibrium (i.e., their basic parameters R_e, L_Δλ, M_s, (M_s/L_Δλ), σ, and I_e, Δλ are such that the strict virial condition is verified). The situation is likely transient because both internal and/or external events may alter one or more parameters so that the strict virial condition is no longer verified. Since the recovery time can vary greatly from galaxy to galaxy, this ideal situation has an ample range of occurrence probabilities, from frequent in some galaxies to never in others. The parameter β can be taken as the signature of how far the system is from full virialization.

The next step will be that of testing the proposed framework at much high redshift (D’Onofrio & Chiosi 2023).

Acknowledgments

The authors thank the anonymous referee for his/her suggestions and comments.

References

All Tables

All Figures

Fig. 1. The size-mass and Faber–Jackson relationships. Left panel: stellar mass vs. radius relations for the Illustris-1 (open red squares) and the Illustris-TNG-100 (blue dots) samples at z = 0 and comparison of the models with the WINGS data (black dots). There are 2400 objects for the Illustris-1 sample and about 600 objects for the Illustris-TNG-100. The mean radii of Illustris-1 are smaller by about a factor of two for stellar masses lower than about 6 × 1010 M⊙, while they are nearly equal (if not slightly larger) above this limit. Right panel: L − σ plane of the same data. The symbols and color-coding are the same as in the left panel.

In the text

	Fig. 2. Comparison between observed and calculated parameters. The data from the WINGS database are used. The black solid line gives the 1:1 relationship. The red solid line is the bi-linear least-squares fit of the distribution.
In the text

Fig. 3. The distribution of the β parameter. Upper panel: histogram of the β solutions derived from Eq. (14) for the real ETGs (black line) compared with the distribution derived for the galaxies of the Illustris-1 simulation when β is calculated looking at the variation in luminosity and velocity dispersion in two close redshift epochs (z = 0.2 and z = 0). Lower panel: histogram of the β solutions derived from Eq. (14) for the real ETGs (black line) and the Illustris-1 galaxies (red line) that are close to the virial equilibrium. The solid black line gives the average value of β.

In the text

	Fig. 4. The β parameter as proxy of virialization. The plot shows how β changes as a function of the degree of virialization (see Eq. (14)).
In the text

	Fig. 5. Histogram of the β solutions derived from Eq. (14) for the real ETGs (black line) of WINGS, the Illustris-1 galaxies (red line) that are close to the virial equilibrium, and the Illustris-TNG galaxies (blue line) in the same condition. The three histograms are qualitatively similar, thus confirming that our analysis depends little on the choice between the two theoretical databases.
In the text

	Fig. 6. Ie − Re plane. The black and red arrows show the direction of motion of galaxies in this plane for large negative and positive values of β. The black solid line gives the least-squares fit of the data, while the broken line represents the zone of exclusion.
In the text

Fig. 7. Ie − Re plane for the WINGS and Illustris-1 data. The black and red arrows show the direction of motion of the WINGS galaxies in this plane for large negative and positive values of β. The green and blue arrows are those of Illustris-1 for negative and positive values of β. In the plot only 1/10 of the Illustris-1 galaxies are used in order to distinguish the objects with different β values moving in different directions.

In the text

	Fig. 8. Re − σ plane. Symbols and colors as in Fig. 6.
In the text

	Fig. 9. Re − σ plane for WINGS and Illustris-1. Symbols and colors as in Fig. 7.
In the text

	Fig. 10. Ie − σ plane. Symbols and colors as in Fig. 6.
In the text

	Fig. 11. Ie − σ plane for WINGS and Illustris-1. Symbols and colors as in Fig. 7.
In the text

	Fig. 12. The log(Re)−log(Ms) plane for the WINGS galaxies. Symbols and colors as in Fig. 6.
In the text

	Fig. 13. The log(Re)−log(Ms) plane for WINGS and Illustris-1 galaxies. Symbols and colors as in Fig. 7.
In the text