© The Authors 2023

1. Introduction

The structural scaling relations (SSRs) of galaxies, that is, the mutual correlations between the main measured structural parameters, such as the effective radius R_e, the effective surface intensity I_e, the total stellar mass M_s, the luminosity L, and the central velocity dispersion σ, were recognized long ago as important tools for understanding the evolution of these stellar systems and for deriving fundamental cosmological information (see e.g., Faber & Jackson 1976; Kormendy 1977; Dressler et al. 1987; Djorgovski & Davis 1987). In particular, the SSRs of early-type galaxies (ETGs), the obtention of which is relatively straightforward, have been used in the past as distance indicators for measuring the Hubble constant (see e.g., Dressler 1987), for testing the expansion of the Universe (see e.g., Pahre et al. 1996), for mapping the velocity fields of galaxies (see e.g., Dressler & Faber 1990; Courteau et al. 1993), and for measuring the variation of the mass-to-light ratio across time (see e.g., Prugniel & Simien 1996; van Dokkum & Franx 1996; Busarello et al. 1998; Franx et al. 2008).

Among the various SSRs, the fundamental plane (FP) relation for the ETGs (Djorgovski & Davis 1987; Dressler et al. 1987), log R_e = a log σ + b log I_e + c, is probably the most studied of the last 30 years. The tilt of this relation with respect to the prediction of virial theorem (VT) has been the subject of many studies (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), invoking different physical mechanisms at work; to mention but a few: (i) the systematic change of 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); (ii) the structural and dynamical non-homology of ETGs (see e.g., Prugniel & Simien 1997; Busarello et al. 1998; Trujillo et al. 2004; D’Onofrio et al. 2008); (iii) the dark matter 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); (iv) the 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); (v) 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); (vi) the effects of dissipationless mergers (Nipoti et al. 2003); (vii) the gas dissipation (Robertson et al. 2006; viii) the irregular sequence of mergers with progressively decreasing mass ratios (Novak 2008); and (ix) the multiple dry mergers of spiral galaxies (Taranu et al. 2015).

A similar long list can be compiled for the small intrinsic scatter of the FP (≈0.05 dex in the V-band), where among the claimed possible physical causes, we have: (1) the variation in the formation epoch; (2) the dark matter content; (3) the metallicity or age trends; (4) the variations of 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; Auger et al. 2010; Magoulas et al. 2012), and so on.

Despite all these efforts, it is still unclear today as to why the FP is so tight and uniform when seen edge-on, while in its projections (i.e., in the I_e–R_e, the I_e–σ, and R_e–σ planes) the distribution of galaxies presents well-defined structures, where regions with large clumps of objects and big scatter are observed together with regions where no galaxies are present, the so-called zone of exclusions (ZOE), and where clear nonlinear distributions are well visible. The mutual dependence of the SSRs, the peculiar shape of the observed distributions, and the link among the various FP projections have never found a single and robust explanation in which the tilt and the scatter of the FP are understood.

The same difficulties are encountered when we consider one particular projection of the FP: the Faber-Jackson (FJ) relation (Faber & Jackson 1976), which is, the correlation observed between the total luminosity L and the central velocity dispersion σ. Even in this case, the observed trend is not that predicted by the VT. In addition to this, it has been shown that the L − σ relation is not consistent with the distribution observed in the I_e − R_e plane, in the sense that it is not possible to transform one space into the other (and vice versa) while adopting the observed classical correlations (see e.g., D’Onofrio & Chiosi 2021).

In other words, there are many underlying questions behind the nature of the SSRs of galaxies. We do not have a single explanation for the tilt of the FP and the shapes of the distributions observed in its projections. We are not able to reconcile the FJ relation and the I_e − R_e plane. We cannot account for the mutual relationship among the various projections of the FP and we do not understand how these planes are mutually linked each other. We do not know wow the SSRs change going back in time and why the FP is so tight.

D’Onofrio et al. (2017) have proposed a new perspective, which may help to simultaneously explain the tilt of the FP and the observed distributions of galaxies in the FP projection planes. The novelty of the approach of these authors is in their assumption that the luminosity of galaxies follows a relation of the form:

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

where t is the time and σ the velocity dispersion. Also, the proportionality coefficient $L_0^{\prime }$ and the exponent β are both functions of time and, more importantly, they can vary from galaxy to galaxy.

This empirical relation is formally equivalent to the FJ relation for ETGs, but 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 the evolution of the stellar content of each object. This new relation mirrors the effects of the evolutionary history of a galaxy on its luminosity and stellar velocity dispersion, parameters that can both vary across time because galaxies evolve, merge, and interact.

In previous papers on this subject, we brought attention to some of the advantages offered by the joint use of the VT and the L = $L_0^{\prime }$(t)σ^β(t) law (D’Onofrio et al. 2017; D’Onofrio et al. 2019; D’Onofrio et al. 2020; D’Onofrio & Chiosi 2021; D’Onofrio & Chiosi 2022; D’Onofrio & Chiosi 2023). Accepting the idea of a variable β parameter, which can take either positive or negative values, yields a simple explanation for the shifts of galaxies along the SSRs. Furthermore, it allows us to understand the physical reasons for the observed distributions of galaxies in the various projection planes. This approach seems to be the correct one because it is able to simultaneously account for: (i) the tilt of the FP, (ii) the existence of the ZOE, and (iii) the shifts of galaxies in the FP projections that are closely connected with the variations of σ and L through the β parameter.

In the present study, we take advantage of what we learned from joining the VT and the law L = $L_0^{\prime }$(t)σ^β(t) to analyze how galaxies move along the SSRs at high redshift. To this aim, and because current observational data at high redshifts are insufficient, we adopt data from the Illustris-1 (Vogelsberger et al. 2014a,b) and the IllustrisTNG (Springel et al. 2018; Nelson et al. 2018; Pillepich et al. 2018a) simulations from z = 0 up to z = 4 and look at the possible changes in the properties of galaxies suggested by the simulations.

The paper is organized as follows: in Sect. 2, we briefly describe the samples of galaxies (both real and simulated) used in our work, we present the basic SSRs at z = 0, and we explain why the simulated data can be trusted at higher redshift. In Sect. 3, we summarize the basic equations of the problem and in Sect. 4, we show how the SSRs change with redshift and how the β parameter is able to account for the observed distributions at each epoch. In Sect. 5, we discuss the β parameter as a thermometer of the virialization condition. In Sect. 6, we discuss the history of mass assembly for a few test galaxies and investigate how β changes as a function of time and history of mass assembly. In Sect. 7, we present our conclusions. Finally, in Appendix A we present a toy model of dry and wet mergers, which we use to estimate the variation of galaxy luminosity as a function of merger and companion star formation.

For the sake of internal consistency with the previous studies of this series, in our calculations with the Illustris-1 database we adopt the same values of the Λ-CDM cosmology used by Vogelsberger et al. (2014a,b): Ω_m = 0.2726, Ω_Λ  =  0.7274, Ω_b = 0.0456, σ₈ = 0.809, n_s = 0.963, H₀  =  70.4 km s⁻¹ Mpc⁻¹. Slightly different cosmological parameters are used for the IllustrisTNG simulations: Ω_m = 0.3089, Ω_Λ = 0.6911, Ω_b = 0.0486, σ₈ = 0.816, n_s = 0.967, H₀  =  67.74 km s⁻¹ Mpc⁻¹ (Springel et al. 2018; Nelson et al. 2018; Pillepich et al. 2018a) . As the systematic differences in M_s, R_e, L, I_e, and σ are either small or nearly irrelevant to the aims of this study, no rescaling of the data is applied.

2. Observational data and model galaxies

The observational data used here are the same as those adopted in our previous works on this subject (see, D’Onofrio & Chiosi 2022; D’Onofrio & Chiosi 2023). The data at redshift z ∼ 0 have been extracted 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 samples used here do not have the same dimensions in each plot because the spectroscopic database is only a subsample of the whole optical photometric sample (containing about 32 700 galaxies). For this reason, in some of our plots we show the distribution of the whole photometric sample, while in others only the subsamples with available measured stellar velocity dispersion or available stellar masses are visible.

The subsample with measured stellar masses M_s contains approximately 1200 galaxies. The masses were estimated by Fritz et al. (2007) by means of the spectral synthesis analysis. This provided the measurements of the stellar masses and of the star formation rate (SFR) at different cosmic epochs (among many other quantities).

The cross-match between the spectroscopic and photometric samples here gives only 480 ETGs with available masses and velocity dispersions. The sample spans a magnitude range from M_V ∼ −16 to M_V ∼ −23 mag, a central velocity dispersion range from σ ∼ 50 to σ ∼ 300 km s⁻¹, and masses from 10^8.5 to 10¹² solar masses¹.

The morphological types of the galaxies were measured with the software MORPHOT for the whole photometric dataset. The final morphological type T is relatively robust, as it comes from the combination of different approaches (see Fasano et al. 2012, for more details).

The error on the measured parameters is ≃20%. These are not shown in our plots because they are much lower than the observed range of variation of the structural parameters in the SSRs or in other words, their size is much smaller than the size of the whole observed distribution of galaxies. Furthermore, no quantitative analysis has been made here, such as fits of data or statistical evaluations.

The sample of real data at z ∼ 0 is used only to demonstrate that the simulated galaxies reproduce the SSRs of the local objects quite well, and therefore there are good reasons to trust the simulation when we look at the behavior of the SSRs at much higher redshift.

The analysis of the SSRs at high redshift is unfortunately still difficult for galaxies above z ∼ 1.0, because the observational surveys at these redshifts contain only a few sparse data points. However, there is some empirical evidence of a varying tilt of the FP with redshift (see e.g., di Serego Alighieri et al. 2005; Beifiori et al. 2017; Oldham et al. 2017; de Graaff et al. 2021).

Given the above difficulties, we decided to perform our analysis of the SSRs at high redshift using the database of artificial galaxies provided by the Illustris-1 and IllustrisTNG simulations. The hydrodynamic simulations, like the Illustris databases, are currently the best available models with which to compare theory with observations, despite the fact that several problems still bias their results.

The first set of artificial galaxies, named Illustris-1, appeared on 2014 (see e.g. Vogelsberger et al. 2014b; Genel et al. 2014; Nelson et al. 2015). Subsequently, a number of works demonstrated that Illustris-1 suffers from a number of problems: it yields an unrealistic population of ETGs with no correct colors, it lacks morphological information, the sizes of the least massive galaxies are too big, and the star formation rates are not always comparable with observations (see e.g., Snyder et al. 2015; Bottrell et al. 2017; Nelson et al. 2018; Rodriguez-Gomez et al. 2019; Huertas-Company et al. 2019; D’Onofrio & Chiosi 2023). In addition to this, there is the claim in the literature that Illustris-1 does 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. 2017; Nelson et al. 2018; Rodriguez-Gomez et al. 2019), while IllustrisTNG produces a much better result (Nelson et al. 2018; Rodriguez-Gomez et al. 2019). There is also the problem of an insufficient number of red galaxies with respect to the observed population of ETGs. For what concerns the internal structure of the Illustris-1 galaxies, Bottrell et al. (2017) measured the Sérsic index, the axis ratio, and the radii, and found that fewer bulge-dominated objects are produced than are seen in observations. In contrast, the IllustrisTNG galaxies have much better internal structural parameters (Rodriguez-Gomez et al. 2019). For this reason, Illustris-1 was superseded in 2018 by IllustrisTNG (Springel et al. 2018; Nelson et al. 2018; Pillepich et al. 2018b).

In this work, we consider only the subsample named IllustrisTNG-100, which is referred to below as simply IllustrisTNG. This sample has approximately the same volume and resolution as Illustris-1 and it used the same initial condition (updated for the different cosmology) adopted by Illustris-1.

Among the many tabulated quantities provided for the galaxies of Illustris-1, we worked in particular with the V-band photometry, the mass, and the half-mass radii of the stellar particles (i.e., integrated stellar populations) for the most massive clusters, for which Cartesian comoving coordinates (x′, y′, z′) are available. In our previous papers, we analyzed the projected light and mass profiles using the z′ = 0 plane as a reference plane. Starting from the V magnitudes and positions of the stellar particles, we computed the effective radius R_e, the radial surface brightness profile in units of r/R_e, the best-fit Sérsic index, and the line-of-sight velocity dispersion. The values of R_e were calculated considering only the star particles inside the friend-of-friends (FoFs) of galaxies and the galaxies inside the FoFs of clusters. We set z′ = 0 to project the coordinates of the stellar particles inside galaxies so that the velocity dispersion is calculated along the z′-axis. The sample does not contain galaxies with masses lower than 10⁹ solar masses at z = 0 because for these objects it was impossible to derive R_e. The total stellar mass has been used here. The dataset for each value of the redshift extracted from the Illustris-1 simulation and used here contains approximately 2400 galaxies of all morphological types. A full description of this dataset is given in Cariddi et al. (2018) and D’Onofrio et al. (2020).

From the TNG-100 dataset, we selected the first 1000 objects – ordered by decreasing stellar mass – provided by the online search galaxy/subhalo catalog². In this case, we used the half-mass stellar radius instead of the effective radius R_e. This radius is not significantly different from the effective radius and its use does not change the conclusions reached here in any way. The data were extracted at redshifts of z = 4, z = 3, z = 2, z = 1, and z = 0 for consistency with those used for Illustris-1.

We chose to use both Illustris-1 and IllustrisTNG for the following reasons: (i) for consistency with our previous works on this subject; (ii) the differences in M_s, R_e, I_e, L, and σ of Illustris-1 and IllustrisTNG do not significantly bias the results on the values of the β and $L_0^{\prime }$ parameters of the L = $L_0^{\prime }$σ^β law (see D’Onofrio & Chiosi 2023; iii) the two data samples are in some way complementary, because IllustrisTNG provides better measurements of the half-mass radii of less massive galaxies, while Illustris-1 contains more massive objects; and (iv) the two simulations agree in terms of the physical parameters of the massive objects.

A detailed analysis of the differences between Illustris-1 and IllustrisTNG data is not addressed here because there are already several studies on this subject (see e.g., Pillepich et al. 2018a,b; Rodriguez-Gomez et al. 2019; Huertas-Company et al. 2019). One major point of tension between the two suites of models concerns the radii of the low-mass galaxies (at roughly M_s ≤ 5 × 10¹⁰ M_⊙, the IllustrisTNG radii are about a factor of two smaller that those of Illustris-1, while above this mass, these radii are approximately equal (Pillepich et al. 2018a,b; Rodriguez-Gomez et al. 2019; Huertas-Company et al. 2019).

Figure 1 shows the distributions of the Illustris-1 (red lines) and Illustris-TNG100 (black lines) data for several parameters used here at two redshift epochs: z = 0 (solid lines) and z = 4 (dashed lines). From the figure, we see that the effective radii of Illustris-1 are systematically slightly larger than those of TNG100. Another significant difference is found in the distribution of the total luminosity and total stellar mass. As mentioned, the Illustris-1 sample does not contain objects with masses of less than 10⁹ solar masses at z = 0. It follows that the distribution of mass versus luminosity appears different for the TNG sample: this distribution is much smoother and flatter than that of Illustris-1, which seems to be peaked at approximately 9 dex for the objects at z = 0. However, the range covered by luminosities and masses is relatively similar. The other parameters appear more or less superposed. We demonstrate below that such differences do not compromise the analysis performed here or our main conclusions. The intrinsic problems of the simulations are of little relevance to our analysis because: (i) we do not make use of the color of galaxies or of the SFRs; (ii) we demonstrate (see, D’Onofrio & Chiosi 2023) that the two samples of Illustris-1 and IllustrisTNG produce very similar distributions of the β and $L_0^{\prime }$ parameters of the L = $L_0^{\prime }$σ^β law; (iii) we show here that the SSRs at high redshift of the two samples are very similar; (iv) the problems of the simulations are mainly related to the inner characteristics of galaxies, such as their shapes, color gradients, etc., that are of little relevance for the point-like view of galaxies adopted here.

Fig. 1. Comparison between the data of Illustris-1 and Illustris-TNG100. The black lines mark the TNG data, while the red ones the Illustris-1 data. The solid line refers to the data at z = 0, while the dashed lines to the data at z = 4. From top left to bottom right we show: the effective radius Re (enclosing half the total luminosity or half the total stellar mass (for TNG); the effective surface brightness Ie; the central star velocity dispersion σ; the stellar mass-to-light ratio (M*/L); the total luminosity (in solar units); and the total stellar mass (in solar units).

For both Illustris-1 and IllustrisTNG, we did not extract information on the morphology of the galaxies. For this reason, ETGs and late-type galaxies (LTGs) are mixed in our comparison plots. This choice originates from the observation that the SSRs of ETGs and LTGs are almost identical. This is clearly seen in Fig. 2, which shows the R_e − M^* (left panel) and the I_e − R_e (right panel) for the ETGs (open circles) and LTGs (filled black circles). The two distributions are very closely superposed in both diagrams. The only exception is that very large R_e are observed only for the most massive ETGs. This is only partially in agreement with Fig. 11 of Huertas-Company et al. (2019), which shows that ETGs and LTGs follow similar trends, but with small systematic differences for the two morphological types. The data of the WINGS database do not suggest any significant difference in the SSRs of LTGs and ETGs. We believe that the effective radii R_e measured in the work of these latter authors are affected by a systematic bias due to the method used to derive R_e. While for the WINGS galaxies, the effective radius was measured as the circle enclosing half the total luminosity, Huertas-Company et al. used the semi-major axis of the best-fitting Sérsic model. This choice likely introduces a systematic effect due to the inclination of the galaxies and the intrinsic shape of the light profiles. In any case, the inclusion of LTGs is a potential source of bias.

Fig. 2. The Re − M* (left panel) and the Ie − Re (right panel) planes for the WINGS galaxy sample. LTGs (black dots) and ETGs (open circles) share the same distributions in these SSRs. The effective radius Re is given in kiloparsecs in the Re − M* plane and in parsecs in the Ie − Re plane. Masses are given in solar units.

In addition, we note that the completeness of the samples is not critical for the conclusions drawn here. Indeed, we do not attempt any statistical analysis of the data, nor do we fit any distribution to derive correlations. The data are only used to qualitatively show how the distribution of galaxies in the various planes can change with the different cosmic epochs and how the L = $L_0^{\prime }$σ^β law and the β parameter can at least qualitatively account for the variations expected or observed across time. The analysis carried out here is indeed somewhat independent of the level of precision reached by the models of the different sources, because we are mainly interested in presenting the method for deciphering the information encrypted in the observational data of the SSRs. The only hypothesis made here is that we can trust the results of simulations at high redshifts. This hypothesis is based on the fact that the simulations are able to reproduce some features of the distributions seen in the FP projections at redshifts z ∼ 0 and the tilt of the FP at z ∼ 1 (see below). The artificial galaxies match the observations relatively well, reproducing the position of the brightest cluster galaxies and the existence of the ZOE. All of the above gives us confidence that the simulations produce galaxies with luminosities, stellar masses, and effective radii that are not significantly different from those of real galaxies.

Figure 3 shows the four main important SSRs for the WINGS and Illustris data. The left upper panel shows stellar mass M^* versus effective radius R_e³. The WINGS data (black dots), the Illustris-1 data (red dots), and the IllustrisTNG data (blue dots) at z = 0 are well visible. We note that the R_e − M^* relation is clearly nonlinear. The galaxies of small mass are distributed almost horizontally, while the brightest galaxies follow a tail with a slope of close to 1. The real and simulated data nicely superimpose each other over the same range of mass, even if the effective radii provided by Illustris-1 for the least luminous galaxies are systematically greater than those observed. In contrast, IllustrisTNG gives much smaller radii for the low-mass galaxies, which we already discussed in the earlier papers of this series (see e.g., D’Onofrio et al. 2020; D’Onofrio & Chiosi 2022; D’Onofrio & Chiosi 2023). Both observations and simulations suggest the presence of the tail for the brightest ETGs, in which radii and masses are almost identical between observations and simulations. The different number of objects in the tail is due to the fact that different volumes are sampled and to the way in which the samples were created: the WINGS and Illustris-1 datasets include only objects from clusters of galaxies, where large ETGs are frequent, while IllustrisTNG takes galaxies from the general field. In addition, the total volume is different in each of the WINGS, Illustris-1 and IllustrisTNG surveys.

Fig. 3. The Re − M* plane and the three different projections of the FP. The black dots are the WINGS observational data. The red and blue dots are the data extracted from Illustris-1 and Illustris TNG-100, respectively. The galaxies with log(Re) < 2.5 are not plotted here in order to better show the bulk of the galaxy distribution.

The upper right panel of Fig. 3 shows the I_e − R_e plane obtained with the same data (here the sample of WINGS galaxies is much smaller than in Fig. 2 because only the subsample is involved). The important takeaway here is that the simulations correctly reproduce the presence of the tail for the brightest ETGs, which is clearly separated from the cloudy distribution of the least luminous galaxies. This tail, already seen in the original paper of Kormendy (1977), has a slope of close to −1 (that predicted by the VT) and has been attributed to the peculiar evolution of the brightest galaxies, which grow in mass by minor mergers (among others, see e.g., Capaccioli et al. 1992). Our conclusion is therefore that both simulations catch the presence of some peculiar features of the I_e − R_e plane: the cloudy distribution of the faint galaxies, the tail formed by the brightest ETGs, and the ZoE, that is, the region totally void of galaxies above the dashed black line in Fig. 2.

The lower panels of Fig. 3 show the other projections of the FP: the I_e − σ and R_e − σ planes. Again, we observe that the data points of the two simulations overlap the observational data relatively closely. In particular, the simulations are able to reproduce the curvature observed in the two distributions.

The good agreement between observations and simulations at z = 0 is a good starting point; it tells us that simulations are able to reproduce the main features of the SSRs at z = 0. However, given that our aim is to use simulations to infer the possible behavior of the SSRs at higher cosmic epochs, we need further proof that simulations are able to catch the structural parameters of galaxies at much higher redshifts. To this end, we used the data of di Serego Alighieri et al. (2005), who analyzed the FP at z ∼ 1. In our Fig. 4, we can appreciate that the FP at this redshift epoch coming from the Illustris data (red and blue circles as before) is in good agreement with the observed FP (black filled circles). The tilt of the plane is identical to a value for the a coefficient a bit lower than 1 (for the Coma cluster, the tilt is a ∼ 1.2). The tilt is different with respect to that measured for the local clusters and this is an indication that there is an evolution of the structural parameters that seems to be reproduced by the simulations. Even in this case, we note that the radii of the galaxies in the simulations are slightly systematically larger than those measured for the real galaxies, but this does not change the FP tilt of the simulated galaxies. Given that the total luminosity of the galaxies is relatively well reproduced by the simulations, the combination of I_e and R_e is likely correct and the different effective radius simply changes I_e in such a way that the galaxy shifts along the FP and not orthogonally to it.

Fig. 4. Fundamental plane at redshift z = 1. The black filled circles are the data of di Serego Alighieri et al. (2005). The red and blue circles are the data from Illutris-1 and IllustrisTNG, respectively. The FP was determined considering only the galaxies with masses in the same range as those of the observational sample (those with M* > 1010 M⊙).

In conclusion, we observe that the artificial galaxies in the simulations are in relatively good agreement with the real galaxies for what concerns the main structural parameters, even at much larger redshifts. The main differences with real galaxies are in terms of stellar content, colors, and SFR, but these differences are not so strong to generate weird distributions in the SSRs. For this reason, we believe that it is possible to extract information on the evolution of galaxies by looking at the distributions of galaxies in the SSRs. When high-redshift data become available in significant quantities, we will be able to better compare observations and simulations and extract useful information on the evolution of galaxies.

3. The basic equations of our framework

Before discussing the main SSRs predicted for furthest cosmic epochs, it is important to summarize the main conclusions drawn by D’Onofrio & Chiosi (2022, 2023) based on the combination of VT and the L = $L_0^{\prime }$σ^β law⁴. This combination is the key novelty of the approach of these authors and understanding the main premise is necessary in order to be able to interpret the results presented in the following sections. The two equations representing VT and the L = $L_0^{\prime }$σ^β law are:

$$\begin{array}{cc}\hfill {\sigma }^2& =\frac{G}{k_v}\frac{M_{\mathrm{s}}}{R_{\mathrm{e}}}\hfill \\ \hfill {\sigma }^{\beta }& =\frac{L}{L_0^{\prime }}=\frac{2\pi I_{\mathrm{e}}{R}_{\mathrm{e}}^2}{L_0^{\prime }}.\hfill \end{array}$$(2)

In these equations, β and $L_0^{\prime }$ are free time-dependent parameters that depend on the particular history of each object. From these equations, we can derive all the mutual relationships existing among the parameters characterizing a galaxy; namely M_s, R_e, L, I_e, and σ. We find:

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

for the I_e − R_e plane, 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 given by

$$\begin{array}{cc}& \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 }}.\hfill \end{array}$$

We then have

$$\begin{array}{c}\hfill 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 }}\end{array}$$(4)

for the I_e − σ relation and

$$\begin{array}{c}\hfill 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 }}\end{array}$$(5)

for the R_e − σ relation. In addition, we have

$$\begin{array}{c}\hfill 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}}\end{array}$$(6)

for the R_e − M^* relation.

It is important to note here that, in all these equations, the slopes of the log relations depend only on β. This means that when a galaxy changes its luminosity L and its velocity dispersion σ, that is, when β has a well-defined value (either positive or negative), the effects of the motion in the L − σ plane are propagated in all the FP projections. In these planes, the galaxies cannot move in all directions, but are forced to move along the directions (slopes) predicted by the β parameter in the above equations. In this sense, the β parameter is the link we are looking for between the FJ (and the FP) and the observed distributions in the FP projections.

In addition, the combination of Eq. (2) gives us another important equation. It is now possible to write an FP-like equation that is valid for each galaxy depending on the β and $L_0^{\prime }$ parameters:

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

where the coefficients

$$\begin{array}{cc}\hfill 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}$$(8)

are written in terms of β and $L_0^{\prime }$. We note that this is the equation of a plane whose slope depends on β and the zero-point on $L_0^{\prime }$. The similarity with the FP equation is clear. The novelty is that the FP is an equation derived from the fit of a distribution of real objects, while here each galaxy independently follows an equation that is formally identical to the classical FP but is of profoundly different physical meaning. In this case, as β and $L_0^{\prime }$ are time dependent, the equation represents the instantaneous plane on which a generic galaxy is located in the log(σ)−log(I_e)−log(R_e) space and consequently in all its projections.

Finally, the combination of the above equations allows us to determine the values of β and $L_0^{\prime }$, the two critical evolutionary parameters. This is possible by writing the following equations:

$$\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}$$(9)

$$\begin{array}{cc}\hfill & \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}$$(10)

Posing now:

$$\begin{array}{cc}& A=log(I_{\mathrm{e}})+log(G/k_v)+log(M_{\mathrm{s}}/L)+log(2\pi )+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 \\ \hfill \end{array}$$(11)

we obtain the following system:

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

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

with solutions:

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

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

The key result here is that the parameters L, M_s, R_e, I_e, and σ of a galaxy fully determine the evolution that is encoded in the parameters β and $L_0^{\prime }$.

Considering the fact that each structural parameter is known with a maximum error of ∼20%, we cannot fully rely on the single values of beta. On average, we instead show that the galaxies move in the SSRs only in the directions defined by beta. Given this premise, we proceed now to show the basic SSRs at much higher redshifts.

4. The SSRs at high redshift

To explore the behavior of the SSRs at high redshifts, we can only relay on simulations because we do not have enough observational data for galaxies at high redshift. Fortunately, despite the small systematic overestimate of the effective radii, the Illustris-1 and IllustrisTNG data are sufficiently precise to be trusted even at high redshifts. Furthermore, as shown by D’Onofrio & Chiosi (2023), both Illustris-1 and IllustrisTNG produce a very similar distribution for the β parameter. Thanks to this, the simulated data can provide reliable insight into the evolution of the SSRs with time. In the following, we show the results for both the Illustris-1 and IllustrisTNG samples currently available to us.

Figures 5 and 6 present the I_e − R_e plane from z = 4 (upper left) to z = 0 (bottom right) for the Illustris-1 and IllustrisTNG, respectively. For Illustris-1, the whole sequence of redshifts is z = 4 (left upper panel), z = 3, z = 2.2, z = 1.6, z = 1.0, z = 0.6, z = 0.2, and z = 0 (right bottom panel) as indicated. In all the panels, the colored dots indicate galaxies with β > 0 and the black points those with β < 0. Crosses and open squares indicate galaxies with SFR greater than the average ⟨SFR⟩ or lower than the average, respectively, at each redshift epoch. The same color code is also used for the arrows indicating the mean slope of β, which was calculated from Eq. (15) for each object of the simulation. This mean value provides approximately the direction of motion in this plane for most of the galaxies, as expected from Eq. (3).

Fig. 5. The Ie − Re plane at different redshift for Illutris-1. From top left to bottom right, we can see the distribution of galaxies at z = 4, z = 3, z = 2.2, z = 1.6, z = 1.0, z = 0.6, z = 0.2, and z = 0. The colored dots mark the galaxies with β > 0 at each epoch, whereas the black symbols mark those with β ≤ 0. The crosses mark the galaxies with SFR greater than the mean ⟨SFR⟩ at each epoch, while the open squares show those with SFR lower than ⟨SFR⟩. The arrows indicate the average direction of motion predicted on the basis of the values of β from Eq. (3); see also Table 1. The colored arrows mark the average for the galaxies with β > 0 and the black arrows those with β ≤ 0. Their length is not related or proportional to any other quantity; it was chosen for graphical reasons.

Fig. 6. The Ie − Re plane at different redshift for the TNG-100 data. From top left to bottom right, we show the distribution of galaxies at z = 4, z = 3, z = 2.2, z = 1.0, and z = 0. The colored dots mark the galaxies with β > 0 at each epoch, whereas the black symbols mark those with β ≤ 0. The crosses mark the galaxies with SFR greater than the mean ⟨SFR⟩ at each epoch, while the open squares show those with SFR lower than ⟨SFR⟩. The arrows indicate the average direction of motion predicted on the basis of the values of β from Eq. (3); see also Table 1. The colored arrows mark the average for the galaxies with β > 0 and the black arrows show those with β ≤ 0. Their length is not related or proportional to any other quantity; it was chosen for graphical reasons.

The sequence of panels indicates that the tails well visible at z = 0 for the brightest galaxies start to appear at z ∼ 1 − 1.5. This epoch probably corresponds to the time in which minor mergers either with or without star formation on already formed massive objects became the typical event, thus increasing both the mass and radius of galaxies (Naab et al. 2009). It is interesting to note that the directions of the arrows, whose slope depends on Eq. (3) (see Table 1), flip progressively with z, in particular for positive β values, assuming the value close to −1 (as predicted by the VT) at approximately z = 0.6, and remain constant thereafter. This slope gives the only possible direction of motion of galaxies in the I_e − R_e plane at each epoch when the evolution proceeds and β changes. The brightest galaxies, which likely reached full virial equilibrium far back in the past, are no longer affected by strong episodes of star formation, and start to move along this direction at z ∼ 1 − 1.5, forming the tail we observe today. Notably, even the galaxies with β ≤ 0 progressively reach the same slope. This happens because several objects have large negative β values. As we show below, both positive and negative values of β are possible and, as demonstrated by D’Onofrio & Chiosi (2021), this is a necessary condition for reproducing the L − σ distribution starting from the I_e − R_e distribution (and vice versa).

A further notable finding is that the galaxies with strong SFR (greater than ⟨SFR⟩) have in general a positive β. When an object with negative β appears on top of the distribution, it has a SFR that is greater than ⟨SFR⟩. Only later on, when the present epoch is approached, we start to see galaxies in the upper part of the cloud with negative β and a lower SFR than ⟨SFR⟩ (the black open squares). These objects might be relatively small compact galaxies where star formation has finished (a possible candidate for this class of galaxies could be M32). Notably, the galaxies with the higher surface brightness have positive β at high redshift and only later is this region of the plot populated by objects with negative β and low SFR.

A very similar behavior is observed in Fig. 6 when we use the IllustrisTNG data. As in Fig. 5, the colored dots indicate the galaxies with positive β and the arrows show the mean value of β. With respect to the previous figure, we note that the colored arrows show fewer changes of direction. We attribute this behavior to the different dimensions of the two samples. For the Illutris-TNG data, the arrows always have a slope of close to −1, a fact that indicates the relatively large values for β. In any case, the trend of populating the upper region of high surface brightness with objects with negative β and low ⟨SFR⟩ is confirmed.

One might now ask wonder what it is that produces the tail for the brightest objects in the I_e − R_e plane, as the arrows at all redshift epochs are directed in approximately the same directions. This can be better understood by looking at the other FP projections and taking into account that the physical mechanisms at work in small and large galaxies can be different.

Figure 7 provides a better demonstration of how the changes of β across time determine the motions in the FP projections (the symbols and color code have the same meaning as those in Fig. 5). In the I_e − σ plane, the curvature formed by the brightest galaxies is much more pronounced. We note again that as β increases, the colored arrows point in the direction of the tail that we see at z ∼ 0. In the figure, we also note that the galaxies with negative β are preferentially at the bottom of the cloud distribution and their number decreases up to z ∼ 0.6 − 1.0. They then increase again in number and tend to crowd the top region of the distribution. As before, the upper region of the distribution with high I_e is populated by objects with high SFR at the most remote epochs, and only when approaching z = 0 do galaxies with negative β and low SFR appear on top of the distribution.

Fig. 7. The Ie − σ plane for the Illustris-1 data. Symbols and colors are as in Fig. 5.

Again, the TNG data (Fig. 8) suggest a similar picture for the I_e − σ plane. With these data, the change of direction due to positive and negative values of β is much more evident and we understand that the tail originates when β increases and the galaxies progressively become much virialized.

Fig. 8. The Ie − σ plane for the IllustrisTNG data. Symbols and colors are as in Fig. 6.

Figures 9 and 10 display the R_e − σ plane with the Illutris-1 and Illutris-TNG data, respectively. Again the slopes of the arrows predicted by Eqs. (3)–(5) are in good agreement with those inferred from the observed distribution of real galaxies and explain the tail formed by the bright galaxies. The same can be said for the R_e − M^* plane (Figs. 11 and 12). In both planes, the tail formed by the brightest galaxies stands out clearly. The slope of the tail in the R_e − M^* plane is very close to 1, as predicted by VT. The bottom right panel (at z = 0) shows in particular that objects with large negative or positive values of β begin to climb the tail as soon as their mass exceeds about 10¹¹ M_⊙.

Fig. 9. The Re − σ plane for Illutris-1. Symbols and colors are as in Fig. 6.

Fig. 10. The Re − σ plane for IllustrisTNG. Symbols and colors are as in Fig. 7.

Fig. 11. The Re − M* plane for Illustris-1. Symbols and colors as in Fig. 5.

Fig. 12. The Re − M* plane for IllustrisTNG. Symbols and colors are as in Fig. 6.

As in the I_e − R_e plane, the TNG data exhibit relatively similar mean values for the slopes predicted at the different redshifts. The slope is close to 1 and this means that the majority of the galaxies are relatively well virialized at z = 0. A further notable finding is that the TNG100 data indicate the presence of relatively massive objects with small radii that are not visible in Illustris-1. These objects might be the class of compact massive galaxies with high I_e also visible in Figs. 6 and 8. Notably, we can see that all these objects have negative β and very low SFR; in other words, they are isolated compact massive galaxies where the SF stopped a long time ago.

The picture we illustrate here using the Illustris-1 and Illutris-TNG data clearly reveals a progressive trend in the galaxies toward the full virial equilibrium, as indicated by the slopes of the arrows when |β|→∞. This condition is reached by the most massive galaxies at approximately z = 1.5 − 1.0. In general, the galaxies with the largest radii have β > 0, both in simulations and observations. This behavior is compatible with the predictions of minor mergers in which galaxies might increase their radius without significantly changing their mass or luminosity (Naab et al. 2009; Genel et al. 2018).

Finally, we note that the FJ relation does not change significantly with redshift. As the redshift decreases from z = 4 to z = 0 through six intermediate steps, we observe that some galaxies have negative β and others have positive β. The fits of the observed distributions reveal that the slope of the FJ relation progressively decreases, passing from nearly 4 (at z = 4) to nearly 2 (at z = 0).

One may legitimately ask why the scatter of the FJ relation does not increase with time, given that there are objects that move almost perpendicularly to the trend indicated by the observed distribution. The same trend is visible with the TNG data (not plotted here). We believe that the scatter cannot increase as a consequence of the merger activity because the maximum possible variation in luminosity that a galaxy might experience does not exceed a factor of two (when a galaxy approximately doubles its mass by merging with a similar object of the same mass and stellar content), which in log units corresponds to a factor of ∼0.3, a very small shift compared with the scale spanned by the data values.

To support this statement, in Appendix A we present a toy model predicting the changes that a galaxy with a mass M₁, age T₁, and luminosity L₁ would undergo in terms of total luminosity as a consequence of a merger with another object of mass M₂, age T₂, and luminosity L₂. The event may or may not be followed by star formation engaging a certain amount of gas with mass M₃.

Using reasonable values for the masses and luminosities of the three components (see Eq. (A.5) in Appendix A), we may expect that the total luminosity first increases and then decreases on a timescale that depends on the amount of matter engaged in the burst of activity. In any case, the luminosity evolution is fast up to a few 10⁸ years after the burst (turnoff mass about 3 M_⊙), slows down up to 10⁹ years (turnoff mass of about 2 M_⊙), and then becomes even slower afterwards. The estimated fading rate of the luminosity is about |Δ(log L/L_⊙)| ≃ 0.015 per gigayear (Gyr) and per unit SSP mass, which must be multiplied by 5.8 to get the real fading rate per Gyr (see the SSP database of Tantalo 2005). Consequently, it is very unlikely to catch a galaxy exactly at the time of maximum luminosity. Equation (A.5) allows us to quickly evaluate the effects of mergers with different combinations of masses and ages of the involved galaxies. However, the examples shown in Appendix A demonstrate that, except for the case of a merger between two objects of comparable mass, in which the luminosity and mass of the resulting objects are double the original ones, mergers among objects of different mass and age, and likely undergoing some star formation during the merger, generate objects that do not keep traces of the merger but simply keep the properties (mass and luminosity) of the most massive component. More details are not of interest here.

The main conclusion of this section of our study can be summarized as follows. The hypothesis that VT and the relation of L = $L_0^{\prime }$σ^β work together to govern the evolution of mass M_s, luminosity L, radius R_e, surface brightness I_e, and velocity dispersion σ leads to a coherent and self-consistent explanation of all the scale relations of galaxies, and also provides a reasonable explanation for the tilt of the FP as demonstrated by D’Onofrio & Chiosi (2022, 2023).

5. The important role of β

In order to better understand the effects played by β, it is necessary to think about the possible variations of R_e and I_e when L and σ vary in the L − σ plane. There are six possible changes of L and σ in this plane: σ either decreases, increases, or remains constant, and the same is true for L. The effective relationship between the two variables depends in turn on β; for example, when β is negative, there is not necessarily a decrease in luminosity, and when β is positive, a decrease in luminosity might also occur (see D’Onofrio & Chiosi 2022; D’Onofrio & Chiosi 2023, for a detailed discussion of this topic). The ambiguity in the direction of evolution can only be solved by looking at the movements of the galaxies in the different SSRs, and in particular by observing the behavior of I_e. When the luminosity of a galaxy changes, both the effective radius R_e and the mean effective surface intensity I_e vary. This happens because R_e is not a true physical radius, similarly to the virial radius for example (which depends only on the total mass), but it is the radius of the circle that encloses half the total luminosity of the galaxy. As galaxies have different stellar populations with different ages and metallicity, it is very unlikely that a change in luminosity would not change the whole shape of the luminosity profile and therefore the value of R_e. If the luminosity decreases passively, in general one 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 to each other and that even a small variation in L might result in ample changes of R_e and I_e.

In summary, as already pointed out, the variations of the parameter β with time are responsible of all the changes observed in the FP projections. This means that the FP problem should be considered from an evolutionary point of view, where time plays an important role and the effects of evolution are visible in all the FP projections. The single SSRs are snapshots of an evolving situation. The L = $L_0^{\prime }$σ^β law is informative of this kind of evolution, and can be used to predict the correct direction of motion of each galaxy in the basic diagnostic planes.

We now show that the parameter β changes with the cosmic epochs and that such variations are in turn related mainly to the change of the mean surface intensity due to the natural variation of the star formation activity with time. We show that β tends to be low when star formation is high and vice versa. However, large scatter is seen at all epochs. Furthermore, we show that β increases considerably if and when the galaxy attains the condition of full virialization, that is, when the two variables M_s and R_e combine in such a way as to yield the measured velocity dispersion (i.e., that measured for the stellar content).

Figure 13 shows the β − log(I_e) plane. The dots of different colors represent the galaxies at different redshifts using the same color code as in Fig. 5. From this plot, it is clear that β increases and log(I_e) on average decreases when the cosmic epoch approaches z = 0 (light gray dots). In the remote epochs (z = 4) and up to z ∼ 1.5, we observe an almost linear dependence of β on I_e, in which β ranges from 0 to ∼20. This condition is an indication that, at these epochs, the galaxies are still far from full virialization. The real data of WINGS (black dots) are very precisely superposed over the simulation data, showing a large spread with large positive and negative values of β.

Fig. 13. The β − log(Ie) plane. The colored dots follow the same color code as in Fig. 5, marking the galaxies at different redshift epochs. The black dots are the data of the WINGS database at z = 0.

This behavior of β is connected with the average SFR at the different cosmic epochs. This is clearly seen in Fig. 14, where we note that, when the SFR is high, the values of β are close to 0–10. The large scatter in β starts to be visible with the gray dots at z ∼ 0.6 − 1.0, the same epoch in which we see the tail of high-luminosity galaxies in the SSRs for the first time. Figure 15 also shows that β takes large positive and negative values preferably in galaxies with masses higher than ∼10¹⁰ M_⊙.

Fig. 14. The β − SFR plane. The colored dots and symbols as in Fig. 13.

Fig. 15. The β − log(Ms) (stellar mass) plane. The color code and symbols are the same as in previous figures.

Finally, Fig. 16 shows β versus the quantity [log(2T)] − log(Ω)], which is a proxy for the virial condition, being the difference between the kinetic and potential energy of the stellar systems. The figure clearly indicates that |β| increases, while β can be either positive or negative when the difference between the two energies approaches zero⁵. We note that β remains very close to zero at high redshift. This means that the galaxies are still far from virial equilibrium. In contrast, at low redshifts, the peak of β falls in the intervals 0 to 20 (z = 1) and 0 to 50 (z = 0) with increasing spread toward both high positive and low negative values.

Fig. 16. The β − [log(2T)−log(Ω)] plane. The color code and symbols are the same as in previous figures.

Figure 17 shows the histogram of the number frequency distribution (N/N_tot) of β in the model galaxies at different redshifts. There is not much difference between the histograms of the Illustris-1 and IllustrisTNG samples. The most remarkable features to note are that, (i) at high redshifts (z ≥ 2), the distribution peaks fall in the interval 4 ≤ β ≤ 0 with a small tail of positive values in the interval 0 ≤ β ≤ 4; (ii) the distribution gradually spreads to higher values of |β| at low redshifts (1 and 0), and both positive and negative values of β are present; and (iii) finally, at low redshifts the peaks are visible both in the positive and negative ranges of β values, and |β| can attain large values.

Fig. 17. Histogram of β for model galaxies of the Illustris-1 (red line) and IllustrisTNG (blue line) samples as a function of redshift.

We observe this large dispersion in β because the term 1 − 2A′/A at the denominator of Eq. (15) becomes very close to zero. Consequently, both log($L_0^{\prime }$) and β diverge. According to the direction from which this zero value is approached, one can have either very large and positive or large and negative values for β. As already discussed, this happens when the system is in conditions of strict virialization. The sign of β depends on the particular history of the 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 operational point of view, we may define “state close to strict virialization” when |β|> 20.

Notably, the Illustris-1 and IllustrisTNG models agree very well with the observational data (black dots in these panels). The inclusion of real dynamics and the hierarchical scenario provides the necessary conditions to clearly demonstrate the action of virialization. The hierarchical scenario by mergers, ablation of stars and gas, harassment, secondary star formation, inflation of dimension by energy injections of various kinds, and so on induces strong variations of the fundamental parameters of a galaxy and hence strong temporary deviations from the virial conditions. However, after this happened, the viral conditions are soon recovered over a reasonable timescale, which can be short or long depending on the amount of mass engaged in the secondary star forming activity and the amount of time elapsed since the star forming event took place (see the burst experiments in Chiosi & Carraro 2002; Tantalo & Chiosi 2004). For this reason, detecting systems on their way back to virial equilibrium is likely a frequent event, which would explain the high dispersion seen on the β–I_e plane. The value of β evaluated for each galaxy can provide a useful hint as to 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 suggested by the large number of galaxies with very small or very high β.

6. The history of mass assembly

The Illustris-1 and IllustrisTNG simulations have made clear that the history of the mass assembly of galaxies is not simple, but is characterized by repeated episodes of mass accretion and mass removal. Figures 18 and 19 show the main SSRs and the β − z plane, respectively. Eight single galaxies of different mass and evolutionary history extracted randomly from the sample are displayed in these plots. These galaxies are taken from our Illustris-1 sample. Each galaxy is indicated by a broken line of different color, while the mass assembly history of each object is represented by a series of dots of the same color. Along each line, there are eight points, one for each value of the redshift from z = 4 to z = 0 according to the list already presented in the previous sections. A very similar figure is obtained with the TNG data and these data are therefore not plotted here.

Fig. 18. The SSRs for eight single galaxies. The different colors mark the single galaxies. The positions in the diagram connected by the colored lines mark the distribution at different redshift epochs.

Fig. 19. The β − z plane for eight single galaxies. Colors are as in Fig. 18.

The β values of these eight galaxies at different redshift are shown in Fig. 19 and are very close to zero at every epoch, with the exception of the blue track. The β values can be entered into Table 1 in order to derive the possible directions of motion at each redshift in each of the planes. For example, the yellow and green objects always have values of beta close to zero (slightly positive). These correspond to slopes of around ∼ − 2 in the I_e − R_e plane, ∼2.5 in the R_e − σ plane, and ∼ − 1.5 in the I_e − σ plane. These are therefore objects of the “big cloud” where galaxies can move in every possible direction. The blue galaxy on the other hand reaches relatively high values of β and it is possible to see that its movements are in the direction close to −1 in the I_e − R_e plane.

It is clear from Fig. 18 that the galaxies do not move in the planes of the SSRs in a continuous and uniform way, but rather they randomly change their position at different epochs. In the same figure, we also note that the galaxies with blue, red, brown, and magenta colors are more massive, more luminous, and with higher σ than the others, even at early epochs (z = 4). More importantly, we emphasize that the same galaxies at epochs closer to the present day than z ∼ 1.5 are able to reach both large positive and negative values of β. Their βs start low and gradually increase, and may reach the state of virial equilibrium. In contrast, the less massive and fainter galaxies have always low βs (see Figs. 18 and 19) close to zero and are located in the region of low σ, I_e, and L. The dwarf galaxies never reach the condition of full virialization.

As already pointed out, the condition of full virial equilibrium can be a transient state in the sense that once reached, it cannot be maintained forever if a galaxy undergoes an event such as a merger, stripping, and/or harassment, which may push it away from this condition. However, virial equilibrium can be recovered again on a suitable timescale, which depends on the relative intensity of the disturbing event. For instance, in the case of a merger between two galaxies of comparable mass, which would most likely be accompanied by intense star formation, the resulting system will not be in virial equilibrium and will take a significant amount of time to reach this condition. On the contrary, the merger between two galaxies of significantly different mass, which is likely to be accompanied by modest star formation, will only slightly depart from the virial condition. If so, we expect that after a certain redshift, only the massive galaxies remain unperturbed by mergers and can move toward the condition of strict virialization, while the low-mass ones are still far from this ideal condition. The few objects displayed in Fig. 19 are typical examples of the above situations.

7. Discussion and conclusions

The aim of this paper is to show that combining VT with the L = $L_0^{\prime }$σ^β relation (which is a proxy for evolution, in which β and $L_0^{\prime }$ vary from galaxy to galaxy and over the course of time) has many positive consequences (see D’Onofrio et al. 2017; D’Onofrio et al. 2019; D’Onofrio et al. 2020; D’Onofrio & Chiosi 2021, for earlier studies along this line of thought). The variation of β and $L_0^{\prime }$ with time traces the path followed by each galaxy in the various SSRs. The combination of the L = $L_0^{\prime }$σ^β law and VT yields I_e − R_e, R_e − σ, I_e − σ, and R_e − M^* relations that nicely reproduce the data and more important strongly 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. Importantly, for each galaxy, this system of equations can be used to determine the values of β and $L_0^{\prime }$. With the aid of these relations, we can determine the instantaneous position and direction of motion of a galaxy on the FP and its projection planes. As these equations are limited to ETGs, we name them fundamental equations of galaxy structure and evolution (D’Onofrio & Chiosi 2022; D’Onofrio & Chiosi 2023).

With this study, we show that the Illustris-1 and IllustrisTNG databases provide the basic parameters of galaxies, and that these are in satisfactory agreement with those derived from observational data for the galaxies at z ≃ 0. The given structural parameters indeed reproduce some distinct features observed in the FP projections, such as the tail of the bright ETGs, the ZOE, and the clumps of low-mass objects.

Based on these simulated data, we look at the SSRs at different epochs (from redshift z = 0 up to redshift z = 4), and highlight their expected behavior. In summary, we show that:

1. The SSRs change with time;

2. The variations of the SSRs can be explained by the variation of the β parameter driving the L = $L_0^{\prime }$(t)σ^β(t) law;

3. When β varies with time, a galaxy can move in the SSRs, but only in well-defined directions, which ultimately depend on β. These directions change with time and, toward z = 0, progressively acquire the slope exhibited by the most massive galaxies which lay along the tail of the bright ETGs at z = 0;

4. The parameter β can take both positive and negative values over time. We therefore suggest that the parameter β can be considered as a thermometer that gauges the virialization conditions. Moreover, β can be either large and positive or large and negative when the galaxies are close to virial equilibrium;

5. The only galaxies that can reach the virial state are those that became sufficiently massive (above 10⁹ − 10¹⁰ M_⊙) at high redshifts (z = 4). These are no longer disturbed by the merging events, which become rare events after z ∼ 1.5 and/or in any case are not influential in terms of changing the mass ratio between donor and accretor;

6. Finally, the L = $L_0^{\prime }$σ^β relation can be considered as an empirical way of catching the temporal evolution of galaxies. The values of β (and $L_0^{\prime }$) mirror the history of mass assembly and luminosity evolution of a galaxy.

The conclusion is that SSRs are full of astrophysical information about galaxy evolution.

Acknowledgments

M.D. thanks the Department of Physics and Astronomy of the Padua University for the financial support.

References

Appendix A: Mergers and bursts of star formation

In order to better quantify the effect of a merger on the integrated light of a galaxy, we present here an elementary model of population synthesis. Let us consider two galaxies with total stellar mass M₁ and M₂ and total luminosity L₁ and L₂ (either bolometric or in some passband). The luminosity is generated by the stars already existing in each galaxy. The two galaxies are supposed to merge. The merger event may or may be not accompanied by star formation induced by the merger itself. Let M₃ be the mass of gas (belonging to one of the galaxies or both) that is eventually turned into newly born stars. For simplicity, we consider this event as a unique single stellar population, SSP, of total mass M₃, generating a total luminosity L₃. If no star formation occurs at the merger M₃ = 0 and L₃ = 0.

The total mass of the system is

$$\begin{array}{c}\hfill M=M_1+M_2+M_3\end{array}$$(A.1)

and the ratios between the mass of each component and the total mass are

$$\begin{array}{c}\hfill {\alpha }_1=\frac{M_1}{M}{\alpha }_2=\frac{M_2}{M}{\alpha }_3=\frac{M_3}{M},\end{array}$$(A.2)

(with no star formation M₃ = 0). Let us suppose that the most massive object of the three is M₁, followed by M₂ and M₃ (with M₂ > M₃). Therefore, we have the sequence α₁ > α₂ > α₃.

The total luminosity is

$$\begin{array}{c}\hfill L=L_1+L_2+L_3\end{array}$$(A.3)

and the corresponding ratios between the luminosity of each component and the total luminosity are

$$\begin{array}{c}\hfill h_1=\frac{L_1}{L}{h}_2=\frac{L_2}{L}{h}_3=\frac{L_3}{L}.\end{array}$$(A.4)

With no star formation, L₃ = 0 and h₃ = 0.

As the luminosity of a galaxy depends not only on the mass but also on the age (the luminosity gets fainter with increasing age) and it may undergo large and fast variations in the presence of star formation, the sequence h₁ > h₂ > h₃ can be easily violated. It may be that h₂ > h₁ and h₃ > h₁ and/or h₃ > h₂. In summary, we define the following identities for the total mass and total luminosity, in which the contribution of each component is revealed:

$$\begin{array}{c}\hfill \begin{array}{cc}& M=({\alpha }_1+{\alpha }_2+{\alpha }_3)M\hfill \\ \hfill & L=(h_1+h_2+h_3)L.\hfill \end{array}\end{array}$$(A.5)

Given these premises, we briefly present the key ingredients of our analysis, namely the relation between the mass and luminosity of single stars and time. This provides clues as to the timescale over which the mass and light of the most massive objects in any generation of stars can vary with age. Each generation of stars forms a SSP with a certain abundance of chemical elements. In a SSP, stars are distributed in mass according to some initial mass function (IMF), in which the number of stars per mass interval typically increases with decreasing stellar mass; that is, many more faint stars of low mass than bright stars of high mass. As stars evolve and die, the total light emitted by a SSP decreases and the SSP becomes fainter and redder with time. Finally, galaxies are made of stars born at different epochs and dying at different times. Therefore, the stellar content of a galaxy can be conceived as a mixture of SSPs of different age, which therefore emit different amounts of light of different colors. The total light is the integral of the light emitted by each SSP weighted by the SFR over the whole history of a galaxy, all as a function of time and chemical composition.

In Fig. A.1, we show the time dependence of the mass (left panel) and luminosity (right panel) of single stars (lifetime and luminosity are taken at the stage of the brightest luminosity attained by the long-lived evolutionary stages) of different masses (in the interval 0.6 to 120 M_⊙). In both panels, the colors indicate the metallicity and the ticks along each curve mark the value of the mass; blue is for low metal content and red for high metal content. Figure A.2 shows the luminosity versus time relationship for SSPs of different chemical composition as indicated (left panel), and model galaxies of different mass as indicated (right panel). The left panel shows the luminosity (in solar units) in the V-passband of the Johnson-Bessell system and of the SSPs with different metallicity (solid lines with different colors where blue is for low metal content and red for high metal content). The black solid line is the best fit of the luminosities for metallicities. In the same panel, and limited to the best fit, we also show the luminosity in the B-passband (dashed line). The right panel shows the V-luminosity in the same photometric system of model galaxies with infall over the entire lifetime of the galaxies. This is the luminosity integrated over the many generations of stars formed in the galaxy at a suitable SFR. Chemical enrichment of the gas from which the stars are formed is taken into account. Stellar models, SSPs, and model galaxies are taken from the Padua library of stellar models, isochrones, and SSPs (see Bertelli et al. 1994; Girardi et al. 1996; Bertelli et al. 2008; Nasi et al. 2008, for all details and ample referencing). The model galaxies were calculated by the authors and described in D’Onofrio & Chiosi (2023). The SSPs in use here are for the Salpeter IMF with a slope of x = 2.35 (in number of stars per mass interval), and the SSP mass and luminosity are named M_SSP and L_SSP with M_SSP = 5.826 M_⊙.

Fig. A.1. Stellar mass in M⊙ vs age in years (left panel). The asterisks along each curve are the star masses in the mass interval 0.6 to 120 M⊙. The color code of each curve indicates the chemical composition. Curves are also labeled in accordance with the mass abundance of metals Z according to the law ΔY/ΔZ = 2.5, where Y is the mass abundance of helium and the initial values are Z0 = 0.0001 and Y0 = 0.23. The values of the metalliciies are Z = 0.0001 (black), Z = 0.0004 (blue), Z = 0.004 (green), Z = 0.008 (magenta), Z = 0.020 (red), and Z = 0.050 (dark red). The mass abundance of hydrogen is given by X = 1 − Y − Z. Same as in the left panel but for luminosity vs age (right panel). The luminosity dips at ages above 1 to 1.5 Gyr correspond to the transition masses MHeF. Along each curve, stars more massive than MHeF undergo a quite He-ignition of the core, while those equal or lighter than MHeF undergo core He-flash.

Notable features of these diagrams are:

(i) The luminosity of single stars may vary by more than two orders of magnitude and the lifetime goes from a few million years (Myr) to more that ten billion years (Gyr) as the mass decreases from 120 to 0.6 M_⊙ with little dependence on the initial chemical composition (at least for our purposes). The interval of variation is much narrower for SSPs and is narrower still for galaxy models. In SSPs, this is simply due to the integration over the mass under some IMF; the more numerous low-mass stars of lower luminosity somewhat quench the luminosity of the brighter less numerous stars. In galaxy models, supposing that each star formation event in the time interval dt can be represented by a SSP of suitable chemical composition emitting the total luminosity l_ssp(t), the total light L is given by

$$\begin{array}{c}\hfill L={\displaystyle {\int }_0^{T_G}}\mathrm{\Psi }(t)l_{\mathit{ssp}}(t)dt,\end{array}$$(A.6)

where Ψ(t) is the current rate of star formation in suitable units.

(ii) The remarkable reduction in the luminosity interval passing from single stars to SSPs and model galaxies is simply due to the integration of the contribution of single stars to the total light of the SSPs of a given age and the integration over time of the contributions from the many generations of SSPs weighted by the SFR.

(iii) The luminosity of SSPs varies by about a factor of 100 with little dependence on the chemical composition passing from young to old SSPs. From these SSPs, we can derive the following best fits:

$$\begin{array}{cc}\hfill log{l}_V& =-0.685logT+6.886\hfill \\ \hfill log{l}_B& =-0.835logT+8.278,\hfill \end{array}$$(A.7)

where l_V and l_B are the luminosities in solar units in the V and B passbands of the Johnson-Bessell photometric system and the age T is in years. These luminosities are for the SSP mass of M_SSP = 5.826 M_⊙. To be applied to a galaxy of stellar mass M_s we must apply the transformations

$$\begin{array}{cc}\hfill log{L}_V& =log{l}_V-log{M}_{\mathit{SSP}}+log{M}_{\mathrm{s}}\hfill \\ \hfill log{L}_B& =log{l}_B-log{M}_{\mathit{SSP}}+log{M}_{\mathrm{s}},\hfill \end{array}$$(A.8)

where M_s is in solar units and is identified with the galaxy mass.

(iv) In the model galaxies, initially the luminosity increases by a factor of ten up to a peak corresponding to the maximum of the SFR at about 1.5 Gyr and then decreases by about the same factor down to the present-day value.

(v) As expected, the total luminosity varies by orders of magnitude among galaxies of different mass.

All this has an immediate effect on the maximum variations of the total luminosity in mergers among galaxies that show different mass and intensity of the companion star formation during the merger event. To illustrate this point, we present a few examples below corresponding to typical situations occurring among real galaxies:

Dry mergers (no additional star formation). With this assumption we have always M₃ = 0 and L₃ = 0.

(a) The case of two identical objects (same mass, same stellar content, and same age) with e.g., M₁ = M₂ = 10¹² M_⊙ and luminosity L₁ = L₂ = 10¹¹ L_⊙. The luminosities are derived from the models in Fig.A.2 (right panel) at the present age. In this case, the mass and luminosity of the resulting object are two times the starting value (α₁ = α₂ = 0.5, h₁ = h₂ = 0.5).

Fig. A.2. Luminosity in L⊙ vs time in years of SSPs (left panel) with the Salpeter IMF (in number of stars per mass interval) with exponent -2.5 and total SSP mass of 5.826 M⊙. The chemical compositions of the SSPs and associated color code are the same as in Fig. A.1. Model galaxies with infall calculated by the authors (right panel). The color code indicates the total baryonic mass of the galaxy in M⊙ whose value is also displayed in the label of each curve. The timescale of gas infall is τ = 1 Gyr. The peak in luminosity corresponds to the maximum SFR. The galaxy models are described by D’Onofrio & Chiosi (2023).

(b) The case of two objects with different mass but the same age (old). Suppose M₁ = 10¹² M_⊙, L₁ = 10¹¹ L_⊙, and M₂ = 10¹¹ M_⊙, L₂ = 10¹⁰ L_⊙. The total mass is M = 1.1 × 10¹² M_⊙ and the total luminosity is L = 1.1 × 10¹¹ L_⊙. Therefore, α₁ = 0.909, h₁ = 0.909 and α₂ = 0.091, h₂ = 0.091. The total luminosity is L = (0.909 + 0.091)L, that is, the light is dominated by the most massive galaxy. The merger has increased the total luminosity by about 10%. This is a sort of lower limit because a dry merger between two galaxies of the same age and differing in mass by more than a factor of ten would, under normal circumstances, be undetectable.

(c) In the case of two galaxies with different mass and different age, and therefore different luminosity, the total luminosity may show the effect of the younger object. Consider the following galaxies, of M₁ = 10¹² M_⊙, L₁ = 10¹¹ L_⊙ (the old object), and M₂ = 10¹¹ M_⊙, L₂ = 5 × 10¹⁰ L_⊙ (the young object, taken near the peak value). After merging, the total mass is M = 1.1 × 10¹² M_⊙ while the total luminosity is L ≃ 1.5 × 10¹¹ L_⊙. Therefore, α₁ = 0.909, h₁ = 0.667, and α₂ = 0.091, h₂ = 0.333. The contribution from the less massive but younger object is important, about half that produced by the more massive but older object. Other combinations of masses and ages, and therefore luminosities, can be easily derived from eq.(A.5). Effects due to differences in the mean chemical composition of the stellar contents can be ignored at a first-order approximation.

Wet mergers (with additional star formation). In this case, in eq.(A.5), M₃ ≠ 0 and L₃ ≠ 0. The major difference with respect to the previous cases is that the stellar activity induced by the merger can be approximated by a single giant SSP with its own mass and luminosity. The total light emitted by M₃ is L₃ = (L_SSP/M_SSP)×M₃. Now the interval spanned by the luminosity at aging SSP is much wider than before, and so the effects can be large, depending on the age. For the merging galaxies, we assume M₁ = 10¹² M_⊙, L₁ = 10¹¹ L_⊙, and M₂ = 10¹¹ M_⊙, L₂ = 10¹⁰ L_⊙. Finally, for the mass of the giant SSP simulating the induced star formation event, we adopt M₃ = 10¹⁰ M_⊙. The total mass is M = 1.11 × 10¹² M_⊙. For the associated luminosity, we adopt three values corresponding to a very young, high-luminosity SSP with L₃ = 1 × 10¹¹ L_⊙; an intermediate age, lower luminosity SSP with L₃ = 1 × 10¹⁰ L_⊙; and an old age, low luminosity SSP with L₃ = 1 × 10⁹ L_⊙. In the first case, the total luminosity is L = 2.1 × 10¹¹ L_⊙, whose relative components are h₁ = 0.476, h₂ = 0.048, and h₃ = 0.476. The burst of star formation contributes to half of the total light and 1% of the total mass. This is a rapidly transient situation that fades on a short timescale. In the second case, the total light is L = 1.2 × 10¹⁰ L_⊙, and the relative contributions to the luminosity are h₁ = 0.834, h₂ = 0.083, h₃ = 0.083; the burst of star formation contributes to about 10% of the light, that is equal to the contribution of the less massive old galaxy. In the last case, the SSP is very old; the three relative contributions are h₁ = 0.908, h₂ = 0.091, and h₃ = 0.009. The occurrence of the burst is nearly undetectable. Other combinations of the parameters can be tested with similar results.

What we learn from these simple tests is that in the case of dry mergers, but for the merger between objects of similar mass in which the mass and light are increased by a factor of two, the effect of mergers among objects with different masses and ages barely affect the light of the originally dominant object. Wet mergers with induced star formation more efficiently leave their fingerprints on the light of the post-merger galaxy. Unfortunately, the bright phase is of short duration —the timescale required by the SSP to evolve from a turnoff in the range of bright massive stars (say 20 M_⊙ with a lifetime of a few Myr) down to a turnoff in the range of low-mass stars (say below 2 M_⊙ with a lifetime of about 1 Gyr); see the panel of Fig. A.1. Therefore, it is by far more probable to catch a galaxy that has undergone a wet merger when the burst of star formation has already faded down to low luminosities. This is an interesting result that could explain why the Faber-Jackson relation and the fundamental plane we see today show little scatter in the observational distribution of galaxies.

All Tables

All Figures

Fig. 1. Comparison between the data of Illustris-1 and Illustris-TNG100. The black lines mark the TNG data, while the red ones the Illustris-1 data. The solid line refers to the data at z = 0, while the dashed lines to the data at z = 4. From top left to bottom right we show: the effective radius Re (enclosing half the total luminosity or half the total stellar mass (for TNG); the effective surface brightness Ie; the central star velocity dispersion σ; the stellar mass-to-light ratio (M*/L); the total luminosity (in solar units); and the total stellar mass (in solar units).

In the text

	Fig. 2. The Re − M* (left panel) and the Ie − Re (right panel) planes for the WINGS galaxy sample. LTGs (black dots) and ETGs (open circles) share the same distributions in these SSRs. The effective radius Re is given in kiloparsecs in the Re − M* plane and in parsecs in the Ie − Re plane. Masses are given in solar units.
In the text

	Fig. 3. The Re − M* plane and the three different projections of the FP. The black dots are the WINGS observational data. The red and blue dots are the data extracted from Illustris-1 and Illustris TNG-100, respectively. The galaxies with log(Re) < 2.5 are not plotted here in order to better show the bulk of the galaxy distribution.
In the text

	Fig. 4. Fundamental plane at redshift z = 1. The black filled circles are the data of di Serego Alighieri et al. (2005). The red and blue circles are the data from Illutris-1 and IllustrisTNG, respectively. The FP was determined considering only the galaxies with masses in the same range as those of the observational sample (those with M* > 1010 M⊙).
In the text

Fig. 5. The Ie − Re plane at different redshift for Illutris-1. From top left to bottom right, we can see the distribution of galaxies at z = 4, z = 3, z = 2.2, z = 1.6, z = 1.0, z = 0.6, z = 0.2, and z = 0. The colored dots mark the galaxies with β > 0 at each epoch, whereas the black symbols mark those with β ≤ 0. The crosses mark the galaxies with SFR greater than the mean ⟨SFR⟩ at each epoch, while the open squares show those with SFR lower than ⟨SFR⟩. The arrows indicate the average direction of motion predicted on the basis of the values of β from Eq. (3); see also Table 1. The colored arrows mark the average for the galaxies with β > 0 and the black arrows those with β ≤ 0. Their length is not related or proportional to any other quantity; it was chosen for graphical reasons.

In the text

Fig. 6. The Ie − Re plane at different redshift for the TNG-100 data. From top left to bottom right, we show the distribution of galaxies at z = 4, z = 3, z = 2.2, z = 1.0, and z = 0. The colored dots mark the galaxies with β > 0 at each epoch, whereas the black symbols mark those with β ≤ 0. The crosses mark the galaxies with SFR greater than the mean ⟨SFR⟩ at each epoch, while the open squares show those with SFR lower than ⟨SFR⟩. The arrows indicate the average direction of motion predicted on the basis of the values of β from Eq. (3); see also Table 1. The colored arrows mark the average for the galaxies with β > 0 and the black arrows show those with β ≤ 0. Their length is not related or proportional to any other quantity; it was chosen for graphical reasons.

In the text

	Fig. 7. The Ie − σ plane for the Illustris-1 data. Symbols and colors are as in Fig. 5.
In the text

	Fig. 8. The Ie − σ plane for the IllustrisTNG data. Symbols and colors are as in Fig. 6.
In the text

	Fig. 9. The Re − σ plane for Illutris-1. Symbols and colors are as in Fig. 6.
In the text

	Fig. 10. The Re − σ plane for IllustrisTNG. Symbols and colors are as in Fig. 7.
In the text

	Fig. 11. The Re − M* plane for Illustris-1. Symbols and colors as in Fig. 5.
In the text

	Fig. 12. The Re − M* plane for IllustrisTNG. Symbols and colors are as in Fig. 6.
In the text

	Fig. 13. The β − log(Ie) plane. The colored dots follow the same color code as in Fig. 5, marking the galaxies at different redshift epochs. The black dots are the data of the WINGS database at z = 0.
In the text

	Fig. 14. The β − SFR plane. The colored dots and symbols as in Fig. 13.
In the text

	Fig. 15. The β − log(Ms) (stellar mass) plane. The color code and symbols are the same as in previous figures.
In the text

	Fig. 16. The β − [log(2T)−log(Ω)] plane. The color code and symbols are the same as in previous figures.
In the text

	Fig. 17. Histogram of β for model galaxies of the Illustris-1 (red line) and IllustrisTNG (blue line) samples as a function of redshift.
In the text

	Fig. 18. The SSRs for eight single galaxies. The different colors mark the single galaxies. The positions in the diagram connected by the colored lines mark the distribution at different redshift epochs.
In the text

	Fig. 19. The β − z plane for eight single galaxies. Colors are as in Fig. 18.
In the text

Fig. A.1. Stellar mass in M⊙ vs age in years (left panel). The asterisks along each curve are the star masses in the mass interval 0.6 to 120 M⊙. The color code of each curve indicates the chemical composition. Curves are also labeled in accordance with the mass abundance of metals Z according to the law ΔY/ΔZ = 2.5, where Y is the mass abundance of helium and the initial values are Z0 = 0.0001 and Y0 = 0.23. The values of the metalliciies are Z = 0.0001 (black), Z = 0.0004 (blue), Z = 0.004 (green), Z = 0.008 (magenta), Z = 0.020 (red), and Z = 0.050 (dark red). The mass abundance of hydrogen is given by X = 1 − Y − Z. Same as in the left panel but for luminosity vs age (right panel). The luminosity dips at ages above 1 to 1.5 Gyr correspond to the transition masses MHeF. Along each curve, stars more massive than MHeF undergo a quite He-ignition of the core, while those equal or lighter than MHeF undergo core He-flash.

In the text

Fig. A.2. Luminosity in L⊙ vs time in years of SSPs (left panel) with the Salpeter IMF (in number of stars per mass interval) with exponent -2.5 and total SSP mass of 5.826 M⊙. The chemical compositions of the SSPs and associated color code are the same as in Fig. A.1. Model galaxies with infall calculated by the authors (right panel). The color code indicates the total baryonic mass of the galaxy in M⊙ whose value is also displayed in the label of each curve. The timescale of gas infall is τ = 1 Gyr. The peak in luminosity corresponds to the maximum SFR. The galaxy models are described by D’Onofrio & Chiosi (2023).

In the text