Galaxies' properties in the Fundamental Plane across time
Mauro D'Onofrio Cesare Chiosi
 


   

   

   

   Department of Physics and Astronomy, University of Padova,
              Vicolo Osservatorio 3, I-35122 Padova

              mauro.donofrio@unipd.it;
         
         
             cesare.chiosi@unipd.it
             
              
             

   



 
  
  
   Using the Illustris-1 and IllustrisTNG-100 simulations we investigate the properties of the Fundamental Plane (FP), that is the correlation between the effective radius , the effective surface intensity  and the central stellar velocity dispersion σ of galaxies, at different cosmic epochs.
  
   Our aim is to study the properties of galaxies in the FP and its projections across time, adopting samples covering different intervals of mass. We would like to demonstrate that the position of a galaxy in the FP space strongly depends on its degree of evolution, that might be represented by the β and L'_0 parameters entering the  law.
  
   Starting from the comparison of the basic relations among the structural parameters of artificial and real galaxies at low redshift, we obtain the fit of the FP and its coefficients at different cosmic epochs for samples of different mass limits.  Then, we analyze the dependence of the galaxy position in the FP space as a function of the β parameter and the star formation rate (SFR).
  
   We find that: 1) the coefficients of the FP change with the mass range of the galaxy sample; 2) the low luminous and less massive galaxies do not share the same FP of the bright massive galaxies; 3) the scatter around the fitted FP is quite small at any epoch and increases when the mass interval increases; 4) the distribution of galaxies in the FP space strongly depends on the β values ( on the degree of virialization and the star formation rate).
  
   The FP is a complex surface that is well approximated by a plane only when galaxies share similar masses and condition of virialization.
   

   Galaxies' properties in the Fundamental Plane across time
    M. D'Onofrio
          1
          C. Chiosi1
          
    Received February 5, 2024; accepted ...
============================================================






§ INTRODUCTION


The Fundamental Plane (FP),  the mutual correlation between the effective radius , the effective surface intensity  and the central stellar velocity dispersion σ of galaxies (log R_e = a logσ + b log I_e + c), has been recognized long ago as an important tool for understanding the evolution of galaxies
<cit.>.
In particular, when the FP was studied using the massive early-type galaxies (ETGs), it was recognized as a good distance indicator <cit.>, and as a useful tool for testing the expansion of the Universe <cit.>, for mapping the velocity fields of galaxies <cit.>, and for measuring the mass-to-light ratio variations across time <cit.>. In general the FP and its projections (I_e vs R_e, I_e vs σ, R_e vs σ)  are diagnostic tools of galaxy evolution.

The main characteristics of the FP of nearby galaxies, essentially all ETGs of large masses (M_s ≥ 10^9 - 10^10 M_⊙), are the tilt of the best-fitted plane with respect to the prediction of the virial theorem  (VT) and its small scatter (≈0.05 dex in the V-band).
The origin of the tilt has been discussed in several studies
<cit.>, invoking different physical mechanisms: i) the systematic change of the stellar mass-to-light ratio (M_s/L) <cit.>; ii) the structural and dynamical non-homology of ETGs <cit.>; iii) the dark matter content and distribution <cit.>;  iv) the star formation history (SFH) and initial mass function (IMF) <cit.>; v) the effects of environment <cit.>; vi) the effects of dissipation-less mergers <cit.>; vii) the gas dissipation <cit.>; viii) the non regular sequence of mergers with progressively decreasing mass ratios <cit.>; ix) the multiple dry mergers of spiral galaxies <cit.>. Many authors have confirmed that the observed tilt disappears when the mass FP is considered,  when the parameters used for the fit are marginally dependent on luminosity <cit.>.

Similarly, the small intrinsic scatter of the plane has never found a clear and definitive explanation. The invoked possible effects at play are: 1) the variation of the formation epoch; 2) the dark matter content; 3) the metallicity/age trends; 4) the variations of the mass-to-light ratio M/L; 5) the mixing of morphological types <cit.>. 

Previous studies have primarily focused on the properties of the FP derived from  massive ETGs at low redshifts. These studies often treated the analysis of the galaxy distribution in the FP-space, considering the projections of the plane , , and  as independent topics for which different explanations were often advanced. Within these FP projections, distinct structural patterns emerge, including areas with pronounced clusters of objects and significant scatter, as well as regions devoid of galaxies (referred to as the 'Zone of Exclusions,' or ZOE), and regions displaying non-linear distribution patterns.

Such a wealth of information has never been used to infer a unitary view of the paths followed by galaxies in the various projection planes of the FP-space during  their evolution[Therein after we call FP-space the space of parameters characterizing a galaxy in which the FP is defined.]. Consequently, a comprehensive and coherent explanation  accounting for the interdependence of structural scaling relations, the peculiar shapes of the observational galaxy distributions in the different planes, and the connections between the various FP projections has remained elusive. This gap in understanding has made it difficult to fully grasp the causes of the tilt and scatter of the FP.

The use of the FP-space as cosmological tool for understanding galaxy evolution requires a global view of the properties of the FP and its projections at different cosmic epochs.
Unfortunately, the exploration of the FP at high redshift and for low-mass galaxies has been notably limited. This limitation arises from the considerable effort required by the  observations of faint, low-mass galaxies at high-redshift. These observations require long lasting campaigns and expensive instrumentation. Existing observational surveys at high redshift typically offer sparse data, primarily focused on the largest galaxies. Nevertheless, some empirical evidence has emerged, suggesting the possible  variation of the tilt of the FP with the redshift <cit.>, as well as deviations of faint galaxies from the FP of their more massive counterparts <cit.>.

Today, thanks to the large efforts in producing reliable cosmological simulations, it is possible to guess the properties of galaxies at different times and follow their evolutionary paths. In this context, the FP and its associated projections become valuable diagnostic diagrams, offering crucial insights into a galaxy's evolutionary state and history.

By exploiting the database of model galaxies of Illustris-1 and IllustrisTNG-100, we aim at infering some useful indication about galaxies evolution looking at the FP and its projections at different cosmic epochs. These simulations offer the best investigation tool currently available to us for a successful comparison between theory and observations, despite the fact that some problems still affect the analysis.

The first version of cosmological simulations named Illustris-1 appeared on 2014 <cit.>. Several works have later shown that there are some problems unsolved by this simulation: it yields an unrealistic population of ETGs with no correct colours, it lacks morphological information, the sizes of the less massive galaxies are too large, and the star formation rates are not always comparable with observations <cit.>. 
Illustris-1 does not seem to produce a realistic red sequence of galaxies due to insufficient quenching of the star formation with too few red galaxies <cit.>, the number of red galaxies is in tension with respect to the observed population of ETGs and, for what concern the internal structure of galaxies, the measured Sérsic index, the axis ratio and the radii, are in marginal agreement with observations <cit.>.

Some years later, in 2018, IllustrisTNG <cit.> was released, a new version of the simulation that seems to produce much better results <cit.>,  in better agreement with observations. In IllustrisTNG galaxies have much better internal structural parameters <cit.>, better colors and radii in closer agreement with available data at low redshift.

In our previous works we convincingly demonstrated that these simulations, although not perfect, do reproduce important features observed in the projections of the FP, and even the tilt of the FP and the small scatter around it <cit.>. 

In the present study, using both libraries of galaxy models, we systematically addressed the FP and its projections at different cosmic times. The study of the FP up to z=2 was addressed by <cit.> using the TNG simulations. We will compare our results with theirs.

However, it is worth clarifying that the main focus of the present study is not to determine the most correct fit of the FP across time, but to investigate the distribution on the FP's at different redshifts of galaxies with different I_e, R_e, L, M_s, σ, and associated parameters  β and  L'_0 defined by <cit.>. Since in each galaxy, I_e, R_e, L, M_s, and σ vary with time (redshift), and so do β and  L'_0.  

The ultimate goal is to show that it is possible to infer some useful indications on the degree of evolution of a galaxy by looking at the observed position of it in the FP-space. The goal can be reached thanks to the new method developed by <cit.>. In brief,  they were inspired by the classical Faber & Jackson relation (FJ; <cit.>) linking the luminosity of L with the stellar velocity dispersion σ of galaxies in the nearby Universe at redshift z≃0. The novelty, however, was to consider that galaxies in the course of their evolution first may change both luminosity and  stellar velocity dispersion  for many reasons (star formation, mass acquisition by mergers or stripping events, natural aging of their stellar populations, and so forth), and second for most of time they are in mechanical equilibrium, that is they obey the VT. Based on these grounds, 
<cit.> first generalized the luminosity-velocity dispersion relation by explicitly including the temporal dependence 


    L(t)  = L'_0(t) σ(t)^β(t).

where t is the time, σ the velocity dispersion, and the proportionality coefficient L'_0 and the exponent β functions of time.  Second, they coupled it to the VT, which is a function of M_s, R_e and σ,  to obtain a system of two equations in the unknowns log L'_0 and β whose coefficients are function of the variables characterizing a galaxy (M_s, R_e, L, σ, and I_e). This system of equations is discussed in detail in Appendix <ref> below.

The new empirical relation (<ref>), although formally equivalent to the FJ relation for ETGs, has a profoundly different physical meaning: β and L'_0 are 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. L'_0 is the most important variable parameter of this equation that empirically mirrors the effects of the evolution while β gives important information on the paths followed by each galaxy in the FP-space in the course of time.
These parameters are found to be good indicators of a galaxy's historical mass accretion, star formation, and evolutionary processes, offering an immediate insight into its current stage of evolution. Indeed the β parameter determines the direction of motion of a galaxy in the FP space.
Adopting this new perspective, it was possible to simultaneously explain the tilt of the FP and the observed distributions of galaxies in the FP projections and, at the same time, to understand the real nature of the FJ relation. 

In a series of papers on this subject, we called attention on some of the advantages offered by the joint use of the VT and the  law <cit.>. Accepting the idea of a variable β parameter, taking either positive or negative values, yields an explanation of the movement of galaxies in the FP-space. This approach allows us to simultaneously account for: i) the tilt of the FP, ii) the existence of the ZoE,  and iii) the flip of galaxies in the FP projections. 

In the present study, taking advantage of what we learned from  galaxies at redshfit z≃ 0, we extended the same analysis of the FP and its projections to galaxies at high redshift. Adopting the data of Illustris-1 <cit.> and IllustrisTNG-100 <cit.> from z=0 up to z=4, we attempted to reconstruct the history of galaxies across time, shading some light on the main mechanisms at work during the evolution.

The paper is organized as follows: In Sec. <ref> we briefly describe the samples of galaxies (both real and simulated) used in our work. In Sec. <ref> we present the FP-space at z=0 and we discuss the results for real and simulated galaxies. In Sec. <ref> we analyze the fits of the FP at high redshift for samples including galaxies of different masses. Here we compare our results with those obtained by <cit.>. In Sec. <ref> we show how the FP-space changes when galaxies with different β values are considered. In Sec. <ref> we discuss the problem of the FP scatter and its relation with the β parameter as a thermometer of the virialization condition. In Sec. <ref> we look at the FP-space distribution of the galaxies as a function of their star formation rate (SFR). In Sec. <ref> we address the problem of the determination of the errors that should be associated to the β parameter. In Sec. <ref> we present our conclusions.  Finally, in Appendix <ref> we summarize the key equations of the new theory of the FP and scale relations,  while in Appendix <ref> we present all the calculations relative to the errors discussed in Sec. <ref>. 
  
For the sake of internal consistency with the previous studies of this series, in our calculations with the Illustris-1 database we adopted the same values of the Λ-CDM cosmology used by  <cit.>:
Ω_m = 0.2726, Ω_Λ= 0.7274, Ω_b = 0.0456, σ_8 = 0.809, n_s = 0.963, H_0 = 70.4 km  s^-1  Mpc^-1. Slightly different cosmological  parameters are used by for the IllustrisTNG-100 simulations: Ω_m = 0.3089, Ω_Λ= 0.6911, Ω_b = 0.0486, σ_8 = 0.816, n_s = 0.967, H_0 = 67.74 km  s^-1  Mpc^-1 <cit.>. Since 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 re-scaling of the data has been applied. 



§ OBSERVATIONAL DATA AND MODEL GALAXIES



. The observational data for the real galaxies are the same adopted in our previous works on this subject <cit.>. The data at redshift z∼ 0 have been extracted from the WINGS and Omega-WINGS databases
<cit.>. This sample contains only ETGs, generally with stellar mass M_s > 10^9   M_⊙.

The combination of the WINGS spectroscopic and photometric samples used here contains ∼ 1690 galaxies with measured ,  and σ. A sub-sample of this, made of 270 galaxies, containing the stellar mass, the star formation rates measured by <cit.> and the morphological types given by MORPHOT <cit.> is also used for some more detailed analyses when necessary. The ETGs sample used here takes also advantage of the data for 24 dwarf galaxies studied by <cit.> including objects of lower masses 
(M_s ∼ 10^7 - 10^9   M_⊙). We explicitly mention the use of this dwarf sample in the work.


   
    
The maximum 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 FP-space. The small size of the errors does not  affect the whole distribution of galaxies. The errors are not considered in our fits of the FP, because the range spanned by the coefficients of the fitted planes, when different samples are used, is greater than the uncertainties on the fit induced by them.


. We adopt the galaxy models of the Illustris-1 and IllustrisTNG-100 spanning ample ranges of masses and redshifts. Some preliminary comments on the two are databases are necessary here. 

The sample of model galaxies in Illustris-1 contains ∼ 2400 objects of all morphological types with masses at z=0 larger than 10^9  M_⊙. We used the data in the V-band photometry, the mass and half-mass radii of the stellar particles (i.e., integrated stellar populations), for which Cartesian comoving coordinates are available. This sample is well discussed in all our previous works <cit.>. For this sample
the values of  are  calculated considering the luminosity growth curves of the galaxies that are members of clusters, following the same technique as in the case of real objects. For the galaxies of Illustris-1 sample, the family-tree of each object has been reconstructed so that it is possible to follow each object along the the cosmic time (redshift).  

From the IllustrisTNG-100 dataset we extracted four samples with 1000 objects at increasing redshift from z=0 to z=4, ordered with decreasing stellar masses. The samples have been obtained using the online Search Galaxy/Subhalo Catalog[see https://www.tng-project.org/data/.]. In this case the half-mass stellar radius is used instead of the effective radius . This radius is not so different from the effective radius and its use does not change  the conclusions reached here at all.  The list of progenitors of each galaxy across the cosmic time is in progress.  

The choice of using both Illustris-1 and IllustrisTNG-100 has the following motivations: i) we want to be consistent 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 bias significantly the results  for  the FP fit  and the values of  β and L'_0 in the  law <cit.>; iii)  the two data samples are in some way complementary: IllustrisTNG-100 has better measurements of the half-mass radii of the less massive galaxies, while Illustris-1 is  richer in massive objects; iv) the two simulations agree in the physical parameters of the massive objects and produce very similar distribution of galaxies in the FP-space, apart from the small differences due to the lower  of the dwarf galaxies in IllustrisTNG-100; v) the use of both samples gives an idea of the degree of uncertainty present in our analysis of the FP-space.

It is important to bear in mind the selection criteria applied to model galaxies when analyzing the FP-space. Specifically, it is noteworthy that the Illustris-1 sample gradually excludes low-mass objects (with masses ≤ 10^9   M_⊙) as we approach redshift z=0. In contrast, the IllustrisTNG-100 sample, which is approximately half the size of the Illustris-1 sample, maintains the mass range of galaxies across all redshifts due to a consistent mass-based selection criteria. Therefore, the theoretical deficiency of low mass objects does not reflect a real  scarcity of these galaxies, but rather a consequence of our galaxy selection criteria.

The mass resolution of the models is about 10^6 M_⊙. Therefore galaxies in the low mass intervals (say up to 10^7  M_⊙ or so)  are described by a handful of mass points, consequently their radius R_e and velocity dispersion σ are poorly determined. Galaxies with masses M_s ≤ 10^7  M_⊙ are excluded from our analsysis. 

The detailed analysis of the differences between Illustris-1 and IllustrisTNG-100 has not been addressed here because it has already been presented in other studies on this subject <cit.>. One of the major issues of tension that may be  relevant for the present study is the radius size of the low mass galaxies (M_s ≤ 5.0× 10^10), where the IllustrisTNG-100 radii are about a factor of 2 smaller that those of Illustris-1  while above it they are nearly equal <cit.>. We will see later that for the aims of the present analysis, the effects of such differences do not introduce significant biases in our interpretation of the FP properties.

In any case, in order to understand the range of the structural parameters that are involved in this work, in Fig. <ref> we present the comparison between the data of Illustris-1 (red lines) and IllustrisTNG-100 (black lines) at the redshifts z=0 (solid lines) and z=4 (dashed lines).
It is clear from the figure that the effective radii of Illustris-1 are systematically larger than those of IllustrisTNG-100. Another significant difference is found in the distribution of the total luminosity and total stellar mass. As already explained, Illustris-1 does not contain objects with mass lower than 10^9  M_⊙ at z=0. It follows that the distributions in mass and luminosity are different for the two samples.
As far as the other parameters are concerned, there are also some differences between the two samples. We will see later that such differences in radii, masses,  luminosities, and velocity dispersions  do not compromise the analysis of the FP  as well as the main conclusions. Despite the small differences in the distribution of galaxies in the parameter space, the global behavior is the same for the two simulations and suggest the same conclusions about the physical effects at play in shaping the FP-space. 

Other points of weakness in the two libraries of model galaxies  are of little relevance for our analysis because: i) we do not make use of the galaxy colors; ii) we have demonstrated <cit.> that  Illustris-1 and IllustrisTNG-100 samples produce very similar distributions of the β and L'_0 parameters of the  law; iii) the FP-space at each redshift for the two samples is very similar, not in the sense that the FP coefficients are identical, but in the general distribution of the galaxies in this space; iv) the point mass view of the galaxies adopted here secures that our analysis is not too much affected by the problems affecting  the inner structure of the model galaxies of the two libraries.
In fact  both for Illustris-1 and IllustrisTNG-100 we did not use information on the  morphology of  galaxies. ETGs and late-type galaxies (LTGs) are mixed together. This choice originates from the observation that the FP projections of ETGs and LTGs are almost identical <cit.>. In addition, we will see that the FP defined by  massive ETGs is only  part of a more general distribution of galaxies in the FP-space. In this sense the presence of LTGs in the samples does not alter none of the conclusions previously reached for the massive ETGs.

The last remark concerns the completeness of the samples. This is not critical for the conclusions drawn here. We will see that the fit of the FP is always performed with a large number of galaxies and that the observed differences can be clearly attributed to the characteristics of the sample under analysis.
This is somehow independent of the level of precision reached by the model galaxies of the two  samples. In fact what we want to emphasize  is not the precision reached in the derivation of the FP, but the general behavior of galaxies in the FP-space as a function of their mass and the dependence of the observed distributions on the β parameter.


 


§ FIT OF THE FP AT Z=0




   
 
 The first step to undertake is to set up the FP and its projections at z=0 for the galaxies of our observational sample and the Illustris-1 and IllustrisTNG-100 galaxies to compare them with other FP's in literature.  
 
 
There are different analytical expressions of the FP, we adopt here the following one
 
 
    log R_e = a log (σ) + b < μ_e > + c

 
where   < μ_e > is the mean effective surface brightness expressed in mag/arcsec^2 and all other symbols have their usual meaning and are expressed in suitable units, see also appendix <ref> and <cit.>. 
[Usually the effective specific intensity I_e is expressed in L_⊙/pc^2. However, limited to the FP, we express I_e as <μ_e>, the mean surface brightness, in mag/arcsec^2. ]
The fit of the WINGS data (mainly massive galaxies with log (M_s/M_⊙) > 10) yields the following coefficients of the FP:  a=0.96, b=0.30, and c=-7.66. 
All our 3D fits are based on the statistical regression of the R program (the free software environment for statistical computing and graphics). 

Figure <ref> shows the FP at z=0 (solid black line) for the WINGS sample to which the 24  dwarf ellipticals of <cit.> are added (black dots and asterisks, respectively) and compare it with the data of Illustris-1 (red dots) and IllustrisTNG-100 (blue dots). The comparison  is made by correlating the R_e of a galaxy (either model or real) with the R_e we would obtain from eqn. (<ref>) representing the FP. The dotted line is the  one-to-one correlation (i.e. perfect coincidence between the two values of R_e for each object).  
The four panels plots the distributions for different ranges of the stellar mass. It is clear from the plots that if the coefficients of the FP  derived from the whole WINGS database  are used also for the fainter, lower mass galaxies, these latter deviate from the main trend. This occurs in particular with the IllustrisTNG-100 dataset. The plume-like features pointing toward directions above the plane are well evident and  quite similar in both samples of model galaxies (there is however a small offsets  due to the different values of  in the two databases). For IllustrisTNG-100 there are also objects that deviate in the opposite direction: these are likely objects with very small  and high effective intensity . The opposite for the upward plume-like features (larger  and lower ). This point will be much more clear when the projections of the FP are shown. These features are nearly absent in the last panel of Fig.<ref> at the bottom right where only galaxies more massive than 10^10  M_⊙ are displayed.

It is also important to note that the deviation observed for the low mass model galaxies is present even in the sample of real galaxies. Here we used the small sample of faint ETGs studied by <cit.> (black asterisks). These objects starts to show the progressive shift from the FP of massive galaxies along the same direction suggested by both libraries of theoretical models, that is toward lower values of  (upward shift). The smaller shift shown by real objects with respect to that seen for models is likely due to the small number of dwarfs (24 only) available in our sample and to selection effects: only the brightest dwarf ETGs with similar  have been detected and studied (i.e., have measured σ).
Their smaller range of  is at the origin of the observed deviations.

In summary, the models suggest that the galaxies with masses lower than ∼ 10^9  M_⊙ do not share the same FP of the more massive ETGs. The galaxies with low  extend toward larger , while those with high , should correspond to those with short . This confirms previous results already obtained during the study of the FP <cit.>.

Furthermore, these plume-like features should not depend on the galaxy morphology; we cannot prove this statement because the information on morphology is missing, but we guess that even for a sample of pure ETGs the trend should be  the same. It is only the range spanned by the specific intensity  the ultimate cause of the plume-like features. Since ETGs and LTGs span approximately the same range of  <cit.>, they should also exhibit the same features.



§ THE FP AT HIGH REDSHIFT



<cit.> demonstrated that Illustris-1 and IllustrisTNG-100 are in quite good agreement with the observational data of the FP-space at z=0. The same is also true for the FP at redshift z=1, at least looking at the data of <cit.> at z∼ 0.8, and those at z=2 according to <cit.>.

We need to clarify that passing from the FP at z=0, for which we have used a sample of galaxies selected according to some criteria, the same galaxies are not  traced back to high redshifts, but galaxies that are present  at high redshift are reselected according to the same criteria. In other words we do not use the family-tree becasuse it is not avalable for both galaxy samples (Illustris-1 and IllustrisTNG-100). Fortunately, this is less of a problem because we are not interested in reconstructing the FP for objects belonging to the same family-tree,  but only in discussing  the statistical properties of the FP at different redshifts. In other words, the main results of our analsysis are not signicantly affected by this type of selection. 

Based on the two libraries of model galaxies to our disposal, we examined the FP at different redshifts, using different  mass intervals of galaxies. In addition to this, we compared the edge-on distribution of galaxies using two representations of the FP: i) at all redshifts the FP is supposed to be the one derived at z=0, that is no evolution of the FP; ii) at each redshift a new FP is derived using galaxies at that redshift. The FP changes with redshift (time). 
Also here, all 3D fits are based on the R method (we have already for the case of z=0). 


   

    
   

 

Figure <ref> shows the edge-on view of the FP at different redshifts for all model galaxies of Illustris-1 and IllustrisTNG-100 with mass greater than log(M_s/M_⊙)>7. 
More precisely, in Fig.<ref> (and all following ones with the same layout) for each object we plot the correlation between the true radius R_e (x-axis) and the R_e we would obtain from the fit of the FP (y-axis). The dotted lines is the one-to-one correlation. The FP on display refers to massive galaxies (log (M/M_⊙) ≥ 10 derived from the WINGs data (it is often referred to as the “reference FP“). 

Two groups of ten panels are displayed. In the top group, the FP is supposed not to change with redshift. The FP determined at  z=0 is supposed to hold at all redshifts (no FP evolution with time). In the bottom group the FP is let change with the redshift (FP evolves with time). In each group,  the top row  of panels shows the Illustris-1 data (red dots) while the bottom row of panels shows the IllustrisTNG-100 data (blue dots). In each panel we also display the redshift z, the total number N of galaxies, the coefficients (a, b, and c) of the FP fit,  and the dispersion Δ of galaxies around it. In each panel, the dotted line is the one-to-one relationship, while the black solid  line  is the FP fit associated to the that panel. The redshift goes from z=4 to z=0 in steps of -1.  Note that the reference FP given in each panel is different for the Illustris-1 and IllustriTNG-100 datasets.   

A careful inspection of all panels reveals that: 
(i) When the FP is kept fixed to the z=0 values obtained with the massive galaxies (log(M_s)>10), a progressive deviations from the one-to-one correlation is observed, in particular when the sample of less massive galaxies is considered. 
(ii) In the case of the Illustris-1 galaxies, the FP of massive galaxies is given by a=1.08, b=0.21, c=-5.87. The dispersion Δ slightly decreases with redshift passing from Δ = 0.08 at z=0 to Δ =0.04 at z=4. For IllustrisTNG-100, the FP is given by  a=1.18, b=0.22, c=-6.38, while the dispersion decreases from Δ = 0.1 at  z=0 to Δ =0.02 at z=4. The agreement of the two determinations of R_e is not satisfactory in particular for the low number of massive IllustrisTNG-100 models at high redshift. 
(iii)  The situation improves when the FP is calculated independently for each galaxy sample at each redshift. The FP coefficients vary with redshift and the dispersion around it significantly smaller with respect to those calculated when the FP is fixed. Notably in this case the galaxies are always close to the fitted planes.
(iv) The analysis clarifies that the FP plane of the galaxy models varies with the redshift. The details depend on the  models in use and  the adopted selection criteria. In general, the agreement between the observed distribution and the FP at varying the redshift is poor if the FP holding at z=0 is supposed to hold at all redshift, while the agreement is good if the FP is left vary with redshift. 

The same analysis was repeated changing the lower mass limit for the galaxies. Now they are those with mass greater than log(M_s/M_⊙)>10. The results are shown in the two groups of panels of Fig. <ref>. All other assumptions and symbols meaning are the same of Fig. <ref>. 

This analysis demonstrated that the coefficients of the fitted FPs vary both with the redshift and with the mass interval in use. As in the case of  Fig. <ref>, when the FP coefficients are fixed, the less luminous and less massive galaxies deviate from the plane of the bright massive objects. On the other hand, when at each redshift the galaxy populations in place at that redshift are used to derive the FP pertinent to the epoch under consideration, a new FP is found with a fairly small scatter.
In addition we note that: i) the scatter around the plane decreases with increasing galaxy mass (it is smaller for the massive galaxies), and ii) the rotation coefficient b of the plane changes by a factor of ∼ 2 when high and low mass objects are separately analysed. The main driver for this change of the rotation coefficient is the effective surface intensity  that varies by several order of magnitudes across the sample.

Another point to keep in mind is that the number of galaxies of different masses varies with redshift. There are two sources of variations. One is real and it depends on the progressive increase of the galaxy mass at increasing cosmic time. The other is an artifact of the selection criteria we adopted to set up the sample of galaxies to examine and also the assumptions made for the mass distribution of the calculated models. All this  may somewhat bias the coefficients of the FP we have determined. The FP coefficients  may somewhat depend on the  size of the galaxy sample <cit.>. However, this effect is only marginally relevant here, because the fits for the FP coefficients are always performed using more than about 30 galaxies.

To somehow cope with this uncertainty affecting the FP coefficients, we performed a series of fits of galaxy distributions in the FP-space by randomly extracting  a fixed number of galaxies from each data sample (each time 200 objects were extracted, and the whole fit procedure was repeated 500 times). In this way we can derive the mean values of the FP coefficients and their scatter, thus providing an estimate of the possible variations induced by the bias
[A problem arises when the number of galaxies is lower than 200. This may occur with the most massive objects. In this case, an acceptable way out is to count massive objects more than once.  Clearly, the final value of the coefficients turns out to be somewhat biased by it.].

Figure <ref> shows the average coefficients of the FP obtained for the Illustris-1 sample (left panel) and for the IllustrisTNG-100 sample (right panel) using the above procedure. For each redshift, we calculated the average values of the coefficients and their scatter around the mean. We also tested the effect introduced by varying the mass interval of the selected galaxies. This is indicated by the color code of the plotted dots.

It is soon evident that the variation with redshift of the FP coefficients is different for the two databases of galaxy models. In Illustris-1, the coefficient a increases from z=0, where a∼ 1.0 at z=0 to a∼ 1.4 at z=2 and remains nearly constant thereafter. This is true for all the mass intervals, with exception of the highest masses (red dots), where a slightly departs from the general trend, most likely because of the small number of galaxies at higher redshift. Conversely, in IllustrisTNG-100 the  coefficient a has smaller variations: it is ∼ 1.2 at z=0 and  a bit smaller at higher z. At z=4, the coefficient a is very different for the different mass intervals, increasing progressively as the mass increases. 
In Illustris-1, the coefficient b  is very similar at all cosmic epochs (b∼ 0.2) with a small scatter around the mean value. This is not the case for IllustrisTNG-100, where the scatter of the coefficient b is much larger, reaching a variation of a factor of ∼ 2  when the low and high mass intervals are used. This difference is likely due to the selection of the galaxies that is different for the two samples. Illustris-1 progressively looses the low mass and faint objects that are characterized by low surface brightness. It follows that in this case the values of b, which strongly depend on , are always very similar.
The  coefficient  c (the zero point of the plane) depends on the coefficient a and shows much larger variations in both datasets. The values of the  coefficients c at each redshift are very different for the different samples with different masses. Notably at z=4, Illustris-1 predicts very similar values of c for the samples, while IllustrisTNG-100 gives different values of c for the different samples with different masses.

Considering that the Illustris-1 sample is biased toward large masses at z=0 and that the galaxy radii are systematically much larger, we are inclined to trust more on the results obtained with the IllustrisTNG-100 data. However, the situation is not fully clear with the present data and it requires a deeper analysis in which the bias introduced by the selection effects is taken into account.

In any case, independently on the reliability of the FP coefficients at each redshift, the basic result obtained here from both databases is that at any epoch massive and dwarf galaxies do not share the same FP. This is well clear in particular for the rotation of the plane that is implicit in the different values of the coefficient b when galaxies of different masses are used.


   

However, it is worth noting that at any redshift there is a mean FP for all mass intervals that is always characterized by a  small scatter. This issue is discussed in great detail in the next sections.



 §.§ Comparison with literature data

We compare here our results with those obtained by <cit.>. In order to frame the comparison in the right context we should remind that: 1) The equation of the FP they adopted is 


    log(R_e[kpc])=a+blog(σ[km/s])+clog(I_e[L_⊙/kpc^2])



   

2) The IllustrisTNG-100 data they consider have stellar masses log(M_s)>9.7. This limit was chosen to take into account the mass resolution of the simulation (1.4×10^6 M_⊙).
3) Their photometric data are in the r-band.

The comparison is shown in the two panels of Fig. <ref> for different values of the redshift, that is  z=0 (left panel) and z=1 (right panel).  In the left panel we show the FP's for the WINGS data (black dots), the Illustris-1 simulations (red dots) and the IllustrisTNG-100 simulations  (blue dots) at z=0. The solid lines are the FP's for the three cases whose coefficientes (a, b, and c) are indicated at the top of the panel according to the adopted color code. Finally, at the bottom of the panel we display the coefficients of the FP of <cit.>.  The same comparison is also made in the right panel of Fig. <ref> for models and observational data at redshift z=1 (same color code). WINGS data are of course missing. The result of this comparison is that 
in the case of  the IllustrisTNG-100 data, the coefficients of the fitted distributions (FP's) are quite similar. The small differences can be easily explained in terms of the different number of galaxies, the different bands (V vs r), and the different strategy used for the fit. We conclude that our analysis confirms the results of <cit.>. 



§ THE SCATTER AROUND THE FP


    
At any redshift, the scatter around the FP is always small when the fit of the galaxy sample is made for galaxies at that redshift and no particular mass interval is selected. In other words, the  classical FP discovered and derived for galaxies at z=0 is not the same at all epochs, but most likely it changes with time <cit.>. Extending and using it as it is at z=0 back in time is not correct. Each epoch has its own FP and the scatter around it is small. In this sense only, the FP is a universal property, whose physical origin is still elusive. 

In Figs. <ref>, <ref>, an <ref>  we have presented the FPs,  the comparisons of these with observational data and model galaxies, and quoted the dispersion Δ of data and models around the FPs without discussing the definition of dispersion we have adopted. Given the  FP  represented the straight line of eq. <ref>  (formally the line R:aX +bY+c), we measure the dispersion of data with respect to this line by the mean value of the distance d_p of all data  points with  coordinates X_p, Y_p from the line R:aX +bY+c. 
For any arbitrary point the distance d_p is:


    d_p =  |a X_p + Y_p + c| /√(a^2 +b^2)


and the mean dispersion < Δ > for all the N data points is 

    < Δ > = 1/N∑ d_p

In our analysis we adopt this definition.

Another way of defining the mean dispersion is to look for any data point at the difference between the observational position O  on the FP and the one indicated by P it would take if eq. <ref> is applied. Therefore we have


    < Δ >  = 1/N∑√((O-P)^2)

This definition is used to compare our results with those in literature.

Finally, in order to highlight the side with respect to the reference FP the deviation occurs, in the discussion below we occasionally  use the simple expression:


    < Δ >  = 1/N∑ (O-P)

which may take both positive and negative values. 

Figure <ref> shows, with different colors for the different mass interval of the galaxies from which the FP is derived, the scatter Δ expected at  different redshifts for models galaxies taken from the Illustris-1 (top panel) and IllustrisTNG-100 (bottom panel) libraries. The two data-sets  yield  slightly different results. In the left panels the  Δ's are calculated according eq. <ref>, while in the right panels according to eq.    <ref> to better show that the FP holding at z=0 cannot be extended to higher redshifts. 


   

If we take IllustrisTNG-100 as the reference sample, because it is less biased in mass, the scatter around the FP is systematically lower for cases with a high mass limit. In any case, the mean value of the scatter is quite small (maximum value Δ≤ 0.1) and similar for both libraries of model galaxies. A possible cause of systematic differences, in this case, could be that the number of galaxies of large mass progressively decreases at increasing  redshift.

The situation changes completely if the FP determined at z=0 is  supposed to hold  at all redshifts and it is used to infer the scatter of galaxies around it.  Figure <ref> (right panel) shows the values of Δ we would expect in such a case. The results are presented both for Illustris-1 (red dots) and IllustrisTNG-100 (blue dots) and limited to two values for the galaxy mass limits, namely log (M_s /M_⊙) ≥ 7 (all galaxies are used, top panel) and log (M_s /M_⊙) ≥ 10 (only the massive ones, bottom panel). The scatter is measured by the difference of the observational log(R_e) and the value calculated using the expression of the FP whose coefficients are derived from the galaxies at z=0 with masses larger than 10^10 M_⊙ (the classical FP). We note that now the scatter is much larger than before, both at changing the redshift and varying the mass limit. Furthermore, in the top panel, the mean scatter is negative, because too many objects would acquire too large R_e if calculated with the expression of FP derived for z=0. 
The negative values for the scatter in this case can be explained by the presence of a large number of objects with mass smaller than 10^10 M_⊙ which have radii smaller  than predicted by the adopted FP of reference (for galaxies more massive than  10^10 M_⊙ at z=0, the classical FP).
In contrast, in the lower panel, when the sample includes only galaxies with mass M>10^10  M_⊙, the scatter becomes positive going toward high redshifts. 

The scatter we have estimated is comparable with the results found by <cit.> and <cit.>. Our values are, however, somewhat smaller than those found those authors. The reason resides in the different definitions of scatter that are used. We adopt the distance method given by relation (<ref>) while they used difference method expressed by relation (<ref>).

In any case, it is worth recalling that the scatter we are talking about cannot be mistaken with the scatter deriving from uncertainties affecting both the theoretical simulations <cit.> and observational results <cit.>. The scatter in question here is caused by the uncorrect use of the FP for massive galaxies at z=0 all over the whole mass range and at all redshifts. 

In summary, the analysis of the FP at different redshifts done with the two libraries of model galaxies suggests that: 1) at each cosmic epoch there is a best fitting plane in the FP-space; 2) this plane is not the same across the whole range of masses, but slightly changes with the mass interval; 3) the scatter around the plane is always quite small (likely a bit larger if the galaxy mass interval is large) when the coefficients of the plane are those derived from galaxies at the redshift under consideration; 4) on the contrary, the scatter can be  significantly large if coefficients of the FP  are those derived from galaxies at z=0, in other words the classical FP is used at any redshift. These conclusions agree with and confirm similar results presented in other recent studies <cit.>.



§ THE FP-SPACE AND THE Β PARAMETER



<cit.> and <cit.> already demonstrated that the behavior of galaxies in the FP-space is tightly related to the value of β  in the  law. More specifically,  the direction of motion of galaxies in this space varies with β, because this parameter encodes the effects of galaxy evolution in terms of merging events, star formation, and many other possible physical processes. Aim of this section is to quantify the above statement that the distribution of galaxies in the FP-space changes  with β. 


.  The equation of the FP for this case is 

    log R_e = 0.92 logσ + 0.28 <μ_e> - 7.15  .


when the dwarf galaxies of  <cit.> are not included.  
Figures <ref>, <ref> and <ref> display the FP-space, both in its edge-on view and its projections, for a sample of galaxies at redshift z=0 with masses exceeding 10^7  M_⊙ and β falling in different intervals. The three figures  are organized in groups of four panels (or physical cases). The group leader is the panel showing the FP on its ordinate (the solid black line)   and is labeled with the value of redshift and range of β of the model galaxies in use (e.g. z=0 and -5 < β <5, first top left panel in Fig.<ref>). The dotted line is the one-to-one relationship for comparison. The other panels surrounding the leader (at the top right, immediate bottom left and bottom right)  show the projections of the FP, that is I_e vs R_e, I_e vs σ, and R_e vs σ in the order.  
In each group of four panels the intervals of  β are different. Finally, in each panel we display the observational data of WINGS (black filled squares), the Illustris-1 models (red filled squares) and the IllustrisTNG-100 models (blue filled squares). 

The careful inspection and mutual comparison of the various groups of panels (physical cases) reveals that position and even the occurrence of galaxies in general and/or particular sub-groups of these on the FP and on its projections find tight correspondence in the values and/or range of values spanned by the β parameter. In other words, we found that some features of the galaxy distribution in the afore mentioned planes come and go according to the underlying values of β. 

To illustrate the point let us look at the case of the faint galaxies. In Fig. <ref>, they are the objects with log(R_e) ≃ 0.5  (in Kpc) and -1 < log I_e < 3 (in L_⊙/pc^2)  in the top right panel  of the first group with -5 <β <5 and corresponds to the plume-like structure above the FP shown in the top left panel of the same group.  They disappear moving to the second group of the same figure and same panel but in which β falls in the range 5 <β <10.  
Conversely, as β becomes larger and more negative, it is well visible only the downward deviation from the FP of the bright galaxies.
The same faint objects deviate from the FP in the first panel of the first group but are absent in the same panel of the second group. 
Furthermore, the galaxies with β in the range  -5 to 5, show both objects along the classical edge-on view of the FP (black solid line), that is the upward deviation typical of the faint galaxies. This deviation is not evident for galaxies with β ranging from 5 to 10, while it is the only feature visible when -5<β<0. As β decreases further and reach large negative values, the IllustrisTNG-100 sample begins to exhibit a downward deviation in the edge-on view of the FP. Remarkably, when β falls below -10, the upward deviation disappears.  


   


   


   
 
Continuing the inspection of the various panels in Fig. <ref> and <ref>, the list of features that come and go depending on the range of β's under consideration gets longer and longer. In general we observe that the most massive galaxies are found along the classical FP quite independently from the value of β. The only difference we can note among all these possible situations is in the number of objects present in each interval of β. Notably for -10< β <-5 only one galaxy is found, while in the interval 
-5<β<5 there are  galaxies along the FP and the classical tail of the  plane. For large positive and negative β's, we always find galaxies along the classic FP, but the distributions in the FP projections may change significantly. The tail in the  plane is better visible for low values of β (-10<β<10).
The obvious conclusion is that at z=0,  galaxies more massive than  10^7  M_⊙  are distributed in a different ways in the FP-space according to their mass and evolutionary conditions that are here indicated by  the β parameter. Variations of β correspond to including or excluding  galaxies of different mass and/or different evolutionary conditions. 

Looking for a plausible physical explanation of this otherwise odd behaviour of galaxies at varying β, we suggest here that it  could be ascribed to the effective surface brightness  and its interplay with  . When  is low and  is large, β gets close to 0. The opposite if  is high and/or  is small or normal.
 Equation <ref> in appendix  <ref> clarifies why β can assume different values, either positive or
negative, either small or large. The reason is that L'_0 depends on the quantity 1-2A'/A (see eq. <ref> in appendix  <ref>), that is  on the ratio between 2log(σ) and the combination of  and  in log units. When 1-2A'/A approaches zero,  L'_0 diverges and consequently β does the same. Values of β close to 0 means that 1-2A'/A is significantly different from 0. The brightest galaxies have in many cases large values of β because 1-2A'/A ∼ 0. Many others, in particular the ETGs along the bright tail of the  relation, have values of β∼ 0. These objects are believed to continuously experience minor dry mergers that increase their radius and decrease their effective surface brightness. Such combination yields 1-2A'/A significantly different from zero.

<cit.> argued that  galaxies with β close to 0 are far from  "full virialization", whereas those with large either positive or negative β's are much closer to this condition. A brief comment on what we really meant with “full virialization”, is mandatory here to avoid any misunderstanding. Indeed, galaxies are always close to the mechanical equilibrium (any momentary deviation from this condition by interaction with another is soon wiped off on a short dynamical timescale), that is for most of time they are very close to virialization. However,  when light is used to derive the galaxies' structural parameters, the main effect of it is that their distribution appears tilted in the FP-space with respect to the virial prediction. Several studies have already demonstrated that, when the half-mass radius and surface-mass density are used as parameters, the tilt disappears <cit.>. This means that light introduces an extra information, different from that provided by the mass.
When we observe galaxies with large values of β (either positive and negative), it means that the combination of , , and M_s/L nearly exactly matches 2log(σ), the condition of virialization written using the light structural parameters, those derived from the luminous component of galaxies. The light parameters  and  do not trace perfectly the mass parameters. The term "full virialization" is used here to say that galaxies are in such peculiar condition, where 2log(σ)∼ A. Many galaxies are far from this condition and have β∼ 0. Indeed their combination of  and  does not match such condition (see Fig. <ref>). 



   


  

Figures <ref> and <ref> provide an explanation of the increase of β across time. Many galaxies of larger dimension and mass can reach the condition of "full virialization" (i.e. 1-2A'/A ∼ 0) and their β increases on both negative and positive sides, depending on the sign of L'_0. The remaining  galaxies of similar dimensions and mass with 1-2A'/A ≠ 0 have β∼ 0. 



   


  

. For the sake of completeness, we briefly discuss the case in which the lower mass limit of galaxies is pushed up to log (M_s/M_⊙) = 10). It is obvious that the number of galaxies to disposal drastically decreases (both observational data and model galaxies). 
The equation of the FP for this case is 

    log R_e = 0.89 logσ + 0.27 <μ_e> - 7.04  ,


not much different from the previous one. Figures <ref> and <ref> illustrate the situation. The layout of the figures and their sub-panels is exactly the same as before. Compared to the reference case there is not much to say, but the overall conclusions we reached are similar although more difficult to defend because of the paucity of data. 


. The situation does not change going to high redshifts. For the sake of brevity we limit the discussion to redshifts z=1 and z=4 (and z=0, for comparison).  The lower mass limit of galaxies is now back to  log (M_s/M_⊙) = 7) while the FP is the one of the WINGS data at z=0.  Contrary to what amply discussed in Sec. <ref>, instead of using the FP equation holding good for the new samples of data (observations and models  with fewer massive galaxies),  we kept also here the same equation for the FP of the reference case at z=0. The uncertainty implicit in this approximation is not relevant in the present qualitative context. Furthermore, to simplify the analysis and presentation of the results we split the data in two groups those with β > 0  and those with β < 0. For the sake of comparison we also consider the case z=0,   Figures <ref>, <ref>, and <ref>  shows the FP and its projections at z=0, z=1, and z=4, respectively.  The layout of the figures is always the same. Galaxies with β > 0 are indicated with filled squares those with β <0  with open squares. The color code has remained the same (black: WINGS; red: Illustris-1, blue: IllustrisTNG-100). 
It is clear that objects with β>0 and β<0 do not share the same distribution.


   



   


   

The distribution of the galaxies changes considerably with respect to that of the WINGS galaxies at z∼ 0. Figure <ref> shows that at z=4 the galaxies with mass greater than 10^7  M_⊙ and β>0 are below the FP, while those with β<0 are above the plane of the WINGS galaxies at z=0.

The β parameter changes also with redshift, but does not seem to much affect the scatter around the FP. At z=0 only a modest scatter around the plane is indicated by the model galaxies of both Illustris-1 and IllustrisTNG-100 databases.

In order to check the relevance of the β parameter as a proxy of the degree of evolution (either in terms of galaxy dynamics that in terms of mass accretion and stellar evolution), we present a  plot showing the stellar mass versus the SFR for the Illustris-1 galaxies at z=0.


    

Figure <ref> shows the Main Sequence (MS) relation. Colors mark different β values. It is apparent from the figure that in different intervals of β we can find only galaxies of certain masses. In such intervals galaxies might have different values of the SFR. Note in particular: 1) the absence of several negative β's along the MS at small masses; 2) the small values of β for few big massive galaxies going toward the quenching state of evolution. These objects have likely experienced in the near past episodes of gas accretion and star formation, well testified by the β parameter.

We conclude that β gives an indication of the stellar mass and the SFR of a galaxy. This information, coupled with the fact that β provides the direction of motion of a galaxy in the FP space, gives this parameter the status of a quantitative proxy of galaxy evolution. 


.  Figure <ref> shows the scatter Δ around the FP against β for objects of different masses at different redshifts. The figure displays two groups of sixteen panels. The group at the left side shows the case of the Illustris-1 model galaxies, while the group at the right side does the same for the IllustrisTNG-100 data. In each panel is indicated the low mass limit of galaxies and the redshifts (z= 0, 1, 2, and 4). The red areas visualize the spread while the horizontal, black, solid lines show the mean scatter.  In both groups, the spread of β increases with redshift. Furthermore, while at z=4 only  small values of β are possible, at z=0 both large positive and negative values of β exist. The mean scatter around the plane does not change very much in all the redshift. In any case, examining the data and results shown in this figure, one should always remind  that the number of massive galaxies becomes smaller going toward high redshifts. 



    

Concluding this section,  we can say that the classical FP is far from being the universal planar distribution valid for all stellar systems. 
The distribution of galaxies in the FP-space changes when objects with different physical conditions and evolution (star formation, mergers, etc.) are taken into account. Each object has its own peculiar position in the FP-space and its projection planes according to its peculiar history of mass accretion,  star formation, and luminosity evolution. 
The effects of these  processes are parameterized by the value of β (and L'_0). 

The old idea of a universal FP is approximately true for the distribution of the large and massive ETGs, objects that are well virialized and may have large positive and negative values of β. The fit of these galaxies in the FP-space does not coincide however with the expectation of the theoretical virial plane, because the fit only gives a sort of average value of the possible values of β <cit.>. The light-measured structural parameters do not match the mass-derived parameters.

The Illustris galaxy models suggest that at each redshift a nearly flat distribution in the FP-space is observed for the galaxies in all mass intervals.  In the present framework this simply reflects the fact that the global structural parameters of galaxies are not so far from those achieved by galaxies when the perfect virial equilibrium is reached. The real situation is that galaxies continuously move in the FP-space because of mass accretion, star formation and luminosity evolution, but the mechanical equilibrium is not very far. The β parameter helps us to have an idea of the effects of such evolution <cit.>.

While the high massive galaxies always follow the classical FP (they are in general very close to the virialization and have large β values, either positive and negative), the smaller galaxies, easily deviate from the ideal virial condition, being their structural parameters strongly affected by feedback effects, mergers and all the other possible physical mechanisms at play during evolution. The most sensible parameter is the effective surface brightness  that easily may induce a strong rotation of the FP via its coefficient b. Dwarfs might have very different  and thus easily change their position in the FP-space.
 


§ THE FP-SPACE AND THE SFR



In this section we  tried to understand whether star formation can influence the distribution of galaxies in the FP-space. The SFR of the real WINGS galaxies was measured by <cit.> and it refers to stellar activity that  occurred in the last 20 Myrs. For Illustris-1 and IllustrisTNG-100 we used the SFRs provided by the libraries of model galaxies.


    
    

    

Figure <ref> shows the FP-space for galaxies of all possible βs (-1000<β<1000) and different SFRs. The SFR is measured in M_⊙/yr. The organization of the figure is similar to those already adopted in previous sections. The low mass limit of galaxies is log (M_s/M_⊙) = 7 and the redshift is z=0. The left panels are for SFRs in the interval 0 to 1 M_⊙/yr, while in the right panels the SFR is between 1 to 5 M_⊙/yr. We note that at z=0 when the SFR increases (SFR > 1 M_⊙/yr), the upward extension characteristic of the small model galaxies disappears and only few objects with high mass are visible in the FP-space. This is a selection effect due to the small number of galaxies with high SFR at z=0. Only the galaxies with large masses have the possibility of reach SFR>1 M_⊙/yr.
Note that few dwarf galaxies with large  and SFR>1 M_⊙/yr are still visible below the plane of the large mass galaxies. Most likely, these objects have experienced recent bursts of star formation or mergers or both that have increased the luminosity and the effective surface brightness.

In Fig. <ref> we restrict the interval of β to -10<β<10.  The figure has the same layout  of Fig. <ref>.  Also here we  note that when the SFR<1 M_⊙/yr  all possible masses (>10^7 M_⊙) are present, while when the SFR>1 M_⊙/yr only galaxies with high mass are visible.  Once more, we note that in general the star forming galaxies are situated along the FP and not far from this. 

Going to higher redshifts this behaviour does not change. At any epoch the galaxies with the lowest SFR are those forming the upward tail above the edge-on FP, while those with highest SFR fall in  general below the plane.

In analogy to what already found and commented for the distribution of galaxies in the FP and its projections at varying the parameter  β, choosing the galaxy sample on the basis of the measured SFR generates different distributions of galaxies in the FP-space.
Our conclusion here is that the FP and its projections are good indicators of the evolutionary conditions of galaxies.

There are two obvious remarks that could be made to the above analysis and conclusion. The first one is why we did not follow galaxy by galaxy of the theoretical samples the merger tree of a candidate galaxy for which the progenitors are known across the cosmic time (redshift). In such a case, all the physical variables L, Re, I_e, M_s, σ, β, L'_0, and SFR  would be known (albeit at discrete redshifts) so that the position on the FP at that redshift would be known. This kind of analysis was already made in 
<cit.>, however limited to a small number of objects (about ten) taken from Illustris-1 sample and limited to the correlation SFR-β. The result was that β remained small at all redshifts for most of the objects, with the exception of two of them in which β got very large (both positive and/or negative) only at low redshift. As already anticipated, the merger three is not yet available for the IllustrisTNG-100 suite of models to our disposal (work is in progress). To cope with this, in the present study, we adopted a different strategy, that is to use all the  objects contained in the samples at each redshift, to select them according to  the mass, SFR, and β, and to plot the filtered objects on the FP at the redshft under consideration. In such a case, 
the simulations show that models with different SFR and β occupy different regions of the FP. In this sense β is a tracer of the mean SFR history in galaxies of a certain mass.  
The SFR history is traced on statistical sense and not object by object.

The second objection, is that the SFR taken from the simulation catalogs is an instantaneous measure of the SFR.  If the SFH is bursty because of the mergers, there is no reason to believe that this value should be related to the structural properties of the galaxies, especially at high redshifts when mergers are likely more frequent. Although this objection sounds correct, it is ill posed and may arise misunderstanding. A merger of two galaxies of given mass may cause star formation in addition to that already taking place in a galaxy due to internal reasons (induced versus spontaneous star formation). The induced star formation can be parametrized by two timescales and a specific intensity.  The first timescale is the duration of a merger (the SF trigger), this time is of the order of the dynamical timescale, or the crossing time scale, or the free-fall time scale. Current estimates set this timescale at about a few 10^8 years. Compared to the typical evolutionary timescales of stellar populations in a galaxy, it can be neglected, The star formation trigger can be considered as an instantaneous event. The second timescale  is the time required to evolve the burst down to the typical age of the background stars already in situ. This can be estimated to be about 1 Gyr (i.e. the lifetime to pass from a turnoff mass of 20 M_⊙  to a turnoff mass of about 1.5 M_⊙). Finally, we have the intensity of the burst that can be measured by the amount of gas mass converted into stars. Since a burst of star formation is ultimately detected and measured by the variation (increase) of the luminosity of the galaxy, it is important to know how the total luminosity of a galaxy resulting from a merger  varies as a function of time. The key quantity to look at is the luminosity per unit mass of a generic event of star formation: this luminosity decreases with increasing time, it may vary by a factor of thousand  over the timescale of about 1 Gyr (see above). The fading rate of a star formation event is very rapid at the beginning and it slows down as it gets older. <cit.> and <cit.> have developed a simple model providing the correct final variation in the luminosity caused by a typical merger. In brief, a merger with induced star formation is thought of as three-players game the two merging galaxies of mass and luminosity M_1 and L_1(t) and M_2 and L_2(t), respectively,  plus  the newly born star with mass and luminosity M_3 and L_3(t). The total mass and luminosity are M=∑ M_j and L(t)=∑_j L_j(t)  where j=1,2,3. Introducing the fractionary masses m_j and luminosities l_j(t), the definition of which is obvious, the effect of a merger with new star formation will depend on m_j and l_j(t). In a typical merger of two objects the  mass ratio is M_1:M_2=10:1 (or higher) while M_3/M_1=0.01 (or smaller). Apart from exceptional circumstances (e.g. very recent merger-burst), the majority of mergers would produce very little effect of the total luminosity emitted by the complex.  The only case in which a sizable effect on L would be present is a merger among objects of nearly equal mass (a rare event) or just an ongoing burst of star formation. 
The Illustris simulations provide data at increasing time intervals in general longer than about 1 Gyr, so that it is extremely unlikely to catch a model galaxy at the peak of the SF burst. 
In conclusion, the SFR provided by the model simulation is the value in a situation far from ongoing mergers and/or bursts.   In any case, both continuous and/or bursting SFR change L, I_e, M, R_e and σ,  and in turn β and L’_0 and therefore the position of the galaxy on the FP.  In other words, the position of a galaxy on the FP can hint the value of the underlying β and SFR. 



§ THE ROBUSTNESS OF THE Β PARAMETER



We established that the β parameter is a reliable indicator of a galaxy's historical star formation and evolutionary trajectory. A value of β near zero suggests that the galaxy is significantly far from full virial equilibrium. This can be attributed to various factors such as AGN or SNe feedback effects, galaxy interactions, mergers, or environmental influences.

On the contrary, both large positive and negative values of β imply that the galaxy is much closer to full virial equilibrium, indicating that evolutionary effects play a less significant role.

This leads us to rise question of the robustness of the β parameter obtained from the solution of the fundamental equations of the VT and the  of appendix <ref> 
<cit.>:
We developed  two distinct methods to estimate the uncertainty associated with β. 

In the first approach, we introduced random errors in each of the parameters in eqns. (<ref>) and (<ref>), ranging from approximately 0% to 20%. The mathematical details on the random perturbation of the input parameters of observed galaxies (WINGS data, 479 objects) are given is the appendix <ref>.  Using this technique we generated 200 simulations with these varying errors, and obtained the number frequency distribution of the β parameter, along with its average value and standard deviation. The results are shown in Fig. <ref> which presents the histogram of the expected percentage change of β due to the introduction of these errors (the solid the black line). It is soon evident that β can fluctuate within approximately a 10% interval in both positive and negative directions, with a mean variation of around -2%.

The blue and green lines in the same Fig. <ref> shows the results of our second method, which assesses potential variations of  β due to errors in the structural parameters. This method relies on the error propagation theory and requires the calculation of the partial derivatives of the formal relationships for β and log L'_0  with respect to the parameters M_s, R_e, L, I_e and σ recalling that all other parameters L, R_e and I_e are each other related by the definition of I_e ∝ L/Re^2 and M_s, R_e and σ are related by the VT. The methods requires also the calculations of variances and co-variances. All details of the errors evaluation, analytical passages, and  calculations are provided in the appendix <ref>. Two cases are considered for the random error distribution. In the first case we adopt a Gaussian random distribution mimicking the PSF of the photometry (blue solid line in Fig. <ref>). In the second case, we adopt a flat random distribution of the uncertainties affecting the observational parameters (green histogram in Fig. <ref>).  See the Appendix <ref> for all details. The agreement between data and simulations is fairly good. Notably the gaussian and flat cases give approximately the same uncertainties.

Figure <ref> highlights the expected mutual correlations  among these derivatives. The high correlation is a consequence of these derivatives being linear combinations of the various structural parameters. Furthermore, while the absolute values of these derivatives are high, they tend to cancel each other out, resulting in a percentage variation in β that closely follows  the simulation-based results. This secures that the expected percentage frequencies of the values for β are on solid ground. 


    


    
    


§ DISCUSSION AND CONCLUSIONS


In this study we investigated the evolution of the Fundamental Plane and its projections on the I_e-R_e, I_e-σ, and R_e-σ planes across cosmic times (redshifts). The analysis has used the Illustris-1 and IllustrisTNG-100 libraries together with the observational WINGS data of real ETGs for comparison. Finally, the observational and theoretical data L, R_e, M_s, I_e and σ have been used to derive the theoretical parameters that encode the past evolution of a galaxy, namely L'_0 and β that are derived from coupling the VT with the generalised luminosity-σ relation . 

Before passing to conclusions, we call attention on a few points of this study that could be subject to criticism: 

i) In real galaxies, the light emitted by their stellar populations can be attenuated by the dust component of the their gas content. The same may happen to the light emitted by galaxies traveling across the cosmic medium toward the observer. In other words, internal and external dust can attenuate the light emitted by a galaxy thus affecting the luminosity and all other variables that are luminosity dependent (i.e. I_e, R_e, M_s, SFR, etc.) thus affecting the determination of the FP and β in turn. Unfortunately, the current observational and theoretical data do not take dust into account. This indeed is a point of uncertainty that cannot be easily cured at the present time.

ii) The observational data we are using to compare our theoretical predictions are somewhat old and limited to redshift z≃0. In literature, there are observational data exploring the FP at z>0 that could be used in our analysis, see for instance <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>. However, in this paper, first we used the WINGS data for the sake of continuity with previous papers of this series, second a tight comparison with  observational data at all redshifts was somehow beyond our aims as the main focus was to highlight the potential ability of the  β-L'_0 method in tracing the SFR history of galaxies.

iii) To derive β and L'_0 for the Illustris-1 and IllustrisTNG-100 simulations, we kept the virial factor k_v constant; see the equation of appendix <ref>.
This quantity depends on the Sersic index n that is not given in the Illustris-1 and IllustrisTNG-100 databases.  Considering the typical range spanned by the Sersic index  in real galaxies we have taken the mean value of n = 4. The maximum uncertainty introduced by k_v is about 0.5 in log units.  In any case uncertainties of this order on β are not enough to explain the large range of values spanned by β. Furthermore,  the analysis of the uncertainties on β and L'_0  caused by  uncertainties on all other parameters (L, I_e, R_e, etc.) shows that the overall effects on  β is rather small. Therefore, the large dispersion shown by β ought to be real and due to physical causes (I_e is likely the main driver). 

Given these premises, the main results of this study can be summarized as follows: 

   
      
  * At any cosmic time there is a plane best fitting the distribution of galaxies in the FP-space;
      
      
  * The best fitting plane, i.e. the FP,  is not the same for galaxies of different masses: low mass galaxies do not share the same plane of massive galaxies;
      
      
  * The coefficients of the best fitting plane, i.e. the FP,   change with the mass of the considered objects  and  the redshift;
      
      
  * The best fitting plane obtained at each epoch is  the FP holding at that epoch, the analog of the FP plane inferred from objects at z=0; 
      
      
  * The scatter around the FP is approximately the same for both simulations and does not change significantly with the cosmic epoch;
      
      
  * The scatter around the plane of the large massive galaxies is in general smaller than the scatter of the smaller low mass galaxies;
      
      
  * The distribution of galaxies in the FP-space changes with the mass of the galaxies considered and with the values of β, the parameter that gives a quantitative idea of the evolutionary stage and physical condition reached by a galaxy.
      
      
  * There is a tight correlation between the rate of star formation and the existence and position of galaxies in the edge-on FP. When the star formation is low (roughly lower than 1  M_⊙/yr) galaxies of any mass (above log (M_s/M_⊙) < 7) are visible in the FP-space, while for  SFR > 1  M_⊙/yr only the massive galaxies are visible along the edge-on FP. In any case, the scatter around the plane is small.
      
      
  *  About 50% of galaxies in the sample have the β parameter confined in the interval -20 < β < 20, all the rest falls outside this range (on both positive and negative side). Since galaxies are supposed to be in perfect virial equilibrium (i.e. their I_e, R_e, L, M_S are compatible with the velocity dispersion), they should have very large values of ± |β|. Our results suggest instead that this fraction of objects is still far from such condition. Since the classical mechanical equilibrium expressed by the VT is likely secured for most of time, this implies that the variations in luminosity (and hence I_e) that follow recent episodes of star formation (in which M_s and likely R_e have also changed) are still on the way back to recover the equilibrium condition. In other words the velocity dispersion derived from the combination of the photometric parameters is not exactly equal to the value derived from the Virial Theorem.
 
      
  *  Full understanding of the relationship between  L'_0, β, and fundamental physical phenomena taking place in a galaxy, e.g. mergers, star formation both spontaneous and induced, natural evolution of the stellar content, and others is still missing despite the many available hints. In brief we know that: (i) β depends on the galaxy mass and evolutionary status; (ii) β spans an ample range of values going from say -1000 to +1000, preference goes to the interval -20 to 20; (iii) The variation of β among galaxies is mirrored by their position on the FP. In general  the position of galaxies with β < 0 is different from those with β > 0 (this is in general true at all redshifts); (iv) Furthermore, also the relation between the value of β (sign and absolute value) and position on the FP seems to change with the redshift; (v) More precisely, typically  β≃ 0  at high redshift (z=4) nearly all over the range of I_e, there is an important tail toward large positive beta for log I_e > 2 at z=2; there is a growth of the positive tail and onset of  an important tail of large negative values for log I_e > 1.5 at z=1; finally, we note an important tails of negative and   positive values of β all over the range of log I_e from 0 to 3  at z=0. (vi) β is strongly related to the star formation rate; (vii) Last, in a galaxy the value of β expressed by its stellar content, mass, radius, etc. should also   change with time as indicated by simple galaxy models <cit.>; (viii) Finally, the specific intensity I_e seems to drive the value of β. Despite all these hints, clear understanding of the rules determining the value of β  (and L'_0 in turn) is still missing. Work is progress to hit the target.     
      
      


In conclusion, large scale cosmological simulations of galaxy formation in hierarchical scheme indicate that galaxies  at any epoch are quite close to the mechanical equilibrium and distribute on a nearly planar surface in the FP-space. This surface is well approximated by a plane when the galaxy sample satisfy the condition of “full  virialization”, that in general is achieved by the massive galaxies at the present epoch.

The variation of the FP-space with redshift and with the β parameter gives to these scaling relations the role of cosmological tools. Indeed, the peculiar distributions of galaxies observed at any epoch in this space might now be linked with the physical mechanisms affecting galaxies in their evolution.


      We thank the anonymous referee for his suggestions and remarks.


   aa 
    
 
 



§ THE BASIC EQUATIONS OF THE NEW THEORY



It may be useful to refresh here the formalism and key results reached by   <cit.> and <cit.>. In the following we drop the explicit notation of time dependence for the sake of  simplicity. The two equations representing the VT and the  law are:



    σ^2    =    G/k_vM_s/R_e
    σ^β   =    L/L'_0 = 2π I_e R^2_e/L'_0.



The first equation is the Virial Theorem, where G is the gravitational constant, and k_v a term that gives the degree of structural and dynamical non-homology. The presence of k_v allows us to write M_s (stellar mass) instead of the total mass M_T. k_v is  a function of the Sérsic index n, that is k_v=((73.32/(10.465+(n-0.94)^2))+0.954)) <cit.>. In the two equations, all other symbols have their usual meaning. 
In these equations, β and L'_0 are time-dependent parameters that depend on the peculiar history of each object. 

.
From these equations one can derive all the mutual relationships existing among the parameters M_s, R_e, L, I_e, σ characterizing a galaxy. We find:

(i) The  plane:


    I_e  = Π R_e^γ

where 


    γ=(2/β)-(1/2)/(1/2)-(1/β)


and Π is a factor that depends on k_v, M/L, β, and L'_0. It is given by


    Π  =  [  (2π/L'_0 )^1/β (L/M_s )^(1/2) (k_v/2π G )^(1/2) ]^1/1/2-1/β .



(ii) The  relation: 

 

    I_e =  [ G/k_vL'_0/2πM_s Π^3/γ ]^β-2/1+3/γσ^β-2/1+3/γ  .


(iii) The  relation:


    R_e =  [ G/k_vL'_0/2πM_s/Π ] σ^β-2/3+γ  .


(iv) The  relation:


    R_e =  [ (G/k_v)^β/2L'_0/2π1/Π ]^2(β-2)/β^2-6β+12 M_s^β^2-2β/β^2-6β+12  .


It is important to note  that in all these relationships the slopes  depend only on β. This means that when a galaxy changes its luminosity L,  and velocity dispersion σ, and has a given value of  β (either positive or negative), the effects of this change in the  plane are propagated in all other projections of the FP-space. In these planes the galaxies cannot move in whatever directions, but are forced to move only 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 FP and the observed distributions in the FP projections.

. In addition, the combination of eqs. (<ref>) gives us another important equation. It is now possible to write a FP-like equation valid for each galaxy and depending on the parameters β and L'_0:



    log R_e = a logσ + b <μ>_e + c


where  <μ_e> is the mean surface brightness <I_e> expressed in magnitudes and the coefficients:


    a     =     (2+β)/3 
    
    b     =     0.26 
    
    c     =     -10.0432+0.333*(-log (G/k_v) - log (M/L) 
           -2*log (2π)-log (L'_0))


are written in terms of β and L'_0. We note that this is the equation of a plane whose slope depends on β and the zero-point on L'_0. 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 formally identical to the classical FP, but of profoundly different physical meaning. In this case, since β and L'_0 are time dependent, the equation represents the instantaneous plane on which a generic galaxy is located in the FP-space and consequently in all its projections.

. Finally, the equation system (<ref>) allows us  to determine the values of β and L'_0, the two basic evolutionary parameters. Let us  write the above equations in the following way:


    β [log(I_e)+log(G/k_v)+log(M_s/L)+log(2π)+log(R_e)] + 
    
        + 2log(L'_0) - 2log(2π) - 4log(R_e) = 0 
    βlog(σ) + log(L'_0) + 2log(σ) + log(k_v/G) - log(M_s) + 
    
     - log(2π) - log(I_e) - log(R_e) = 0.



Posing now: 


    A      =    log(I_e)+log(G/k_v)+log(M_s/L)+log(2π)+ 
          log(R_e)  
    
    B      =     - 2log(2π) - 4log(R_e)  
    
    A'     =    log(σ)  
    
    B'     =     2log(σ) - log(G/k_v) - log(M_s) - log(2π) -  
          log(I_e) - log(R_e)

we obtain the following system:


    Aβ + 2log(L'_0) + B = 0 
    
    A'β + log(L'_0) + B'= 0

with solutions:


    β    =    -2log(L'_0) - B/A
    log(L'_0)     =    A'B/A - B'/1-2A'/A.


The key result is that the  parameters L, M_s, ,  and σ of a galaxy fully determine its evolution in FP-space that is encoded in the parameters β and L'_0.
Given this premise,  we proceed now to show the basic scale relationships at  increasing redshift, that is at decreasing time.  



§ THE ERRORS OF Β


 

 §.§ Useful definitions
 
Suppose we have a sample on N_max objects (galaxies), each of which is characterized by the following parameters: luminosity L in some passbands and solar units, stellar mass M_s in solar units, effective radius R_e in kpc or pc, velocity dispersion σ in km/s, and specific intensity I_e = L/2π R_e^2 (where  by convention R_e is expressed in pc). 
 
First of all,  for the sake of convenience, all these parameters are expressed in logarithmic form and in order to simplify the notation they are replaced by  M  = log(M_s), R  = log(R_e), S  = log(σ),  L  = log(L), I  = log(I_e). Finally, a generic object of the sample is indicated by the index i=1,2,......N_max.     
 
With the aid of these parameters, the VT yielding the velocity dispersion as a function of the mass M_s and effective radius R_e, and the relation L = L'_0 σ^β, we set up a system of equations, the solution of which is represented by equation (<ref>) and (<ref>) providing the values of  β and L'_0. It is obvious that β and L'_0 vary from galaxy to galaxy, and for each galaxy with time. Furthermore, since the parameters L, M_s, and so forth,  are surely affected by some observational uncertainty, the natural question arises:  “How  stable are β and L'_0? ”. In addition to it, we considered the characterizing parameters as independent quantities. However this is not the case,  for instance the total luminosity depends on the mass, the determination of the effective radius R_e depends on the luminosity (and mass) profile across the galaxy, the velocity dispersion σ, derived from the virial theorem, depends on the mass and radius, finally the specific intensity I_e depends on the luminosity and the radius. 



 §.§ Simulations of the errors for the galaxies in the sample


The N_max galaxies in the sample have characterizing parameters L, M_s, R_e, σ, I_e whose value has been measured only once and is surely affected by some uncertainty that is difficult to assess for each object. To cope with it and be able to apply the error propagation theory on  β and L'_0 of each object, we simulate artificial errors for each quantity by artificially generating for each galaxies of the list a number of virtual objects whose parameters have been slightly varied with respect to the original ones by some random correction in this way mimicking different sets of measurements. The number of moke measurements is N_err, whose value has to be suitably fixed (say around 50). The range of uncertainty for each value is of the order of 10 to 20% at maximum. Since we are using variables expressed in the logarithmic scale, uncertainties of 10 to 20% correspond to uncertainties of 0.042 to 0.079 in the logarithms, indicated here by  θ. Looking at the mass as an example (we used here the logarithmic variables in compact notation), the procedure is a follows:  given any mass M_0  this is perturbed  by Δ M,  M = M_0 + Δ M so that  


    M = M_0  + Δ M =  M_0 + (ℛ - 0.5) θ


where 
ℛ is a random number from 0 to 1, and θ the adopted maximum uncertainty. In this way the correction Δ M can be positive or negative and the total maximum range for Δ is ±θ M_0.  The same procedure is repeated for R and L. Recalling that the velocity dispersion is derived from the VT and the specific luminosity I_e ∝ L/R_e^2 (the variables are correlated each others), we can simply write


    Δ S = (Δ M - Δ R)/2       and    Δ I = Δ L - 2 Δ R


With new perturbed values we calculate log L'_0 and β. The same procedure  is repeated a number of times (N_err) for each galaxy of the sample and for all galaxies of the list. The resulting frequency (or percentage) distribution of the β's was already shown in Fig.<ref>.  

Another alternative formulation for the errors affecting the parameters is a gaussian distribution mimicking the mathematical behaviour of the point spread function (PSF) used to measure the luminosity profile of the source galaxy out of which the other variables R_e and M_s are derived.
To this aim we adopt  


    F(x) = 1/σ√(2π) exp[- (x-μ)^2/2σ^2]


where F(x) is the normalized distribution relative to the expected ideal case μ and σ the standard deviation (although the same symbol is used, it should not be mistaken with the velocity dispersion), and finally x is the distance both positive or negative with respect to the position of the maximum of F(x). In our case, the obvious choice for μ and σ are μ=0 and σ≃θ.  In order to take the observational uncertainty into account, σ≃ 0.5 or 0.1. With these choices  F(x) ≃ 0 at x=± 0.5 (actually it is zero also at smaller values of x. Any contribution beyond these limits can be ignored for practical purposes.  Therefore, the range of interest is -0.5 ≤ x ≤ 0.5. The maximum F(x) for x=0 is F(0)= 1/σ√(2π).  To determine x, we start from the relation

                          
    x = (ℛ - 0.5)


where ℛ is the random number in the interval 0 to 1,  derive F(x),  use the assigned difference F(0) - F(x) divided by F(0) as the fraction at uncertainty with respect to its maximum value θ, and finally multiply all this by θ to get the real uncertainty at x. Again taking the mass as an example, we alter its value (in logarithmic scale) as follows 


    M = M_0  + Δ M =  M_0 + x/|x|×((F(0) - F(x')/F(0)) ×θ


the uncertainty increases going toward the edges of the interval permitted to x and is zero at the center of the gaussian distribution.  The definition of Δ M is straightforward.  The same procedure is repeated for R and L. Also here, uncertainties for the velocity dispersion sigma Δ S and the specific intensity Δ I are derived from Δ M, Δ R and Δ L. Finally, the whole procedure  is repeated a number of times (N_err) for each galaxy  of the sample and for all galaxies of the list.



 §.§ Elements of the error propagation theory


Let start defining the following auxiliary quantities:

(1) 𝒮_L = ∑_i L_i, 𝒮_M = ∑_i M_i, 𝒮_R = ∑_i R_i, 𝒮_σ = ∑_i σ_i, and 𝒮_I = ∑_i I_i, where i=1,2,......N_err; 


(2) 𝒮𝒮_L = ∑_i (L_i)^2, 𝒮𝒮_M = ∑_i (M_i)^2, 𝒮𝒮_R = ∑_i (R_i)^2, 𝒮𝒮_σ = ∑_i (σ_i)^2, 𝒮𝒮_I = ∑_i (I_i)^2, where i=1,2,......N_err; 


(3) 𝒮𝒮𝒮_ML = ∑_i (M_i× L_i),   𝒮𝒮𝒮_MR = ∑_i (M_i× R_i),
    𝒮𝒮𝒮_MS = ∑_i (M_i× S_i),   𝒮𝒮𝒮_MI = ∑_i (M_i× I_i),
    𝒮𝒮𝒮_LR = ∑_i (L_i× R_i),   𝒮𝒮𝒮_LS = ∑_i (L_i× S_i),
    𝒮𝒮𝒮_LI = ∑_i (L_i× I_i),   𝒮𝒮𝒮_RS = ∑_i (R_i× S_i),
    𝒮𝒮𝒮_RI = ∑_i (R_i× I_i),   𝒮𝒮𝒮_SI = ∑_i (S_i× I_i), where i=1,2,......N_err; 


(4) Averages:   ℳ_L = ∑_i L_i/N_err = 𝒮𝒮_L/N_err, 
                ℳ_M = ∑_i M_i/N_err = 𝒮𝒮_M/N_err,
                ℳ_R = ∑_i R_i/N_err = 𝒮𝒮_R/N_err,  
                ℳ_S = ∑_i S_i/N_err = 𝒮𝒮_S/N_err, and
                ℳ_I = ∑_i M_i/N_err = 𝒮𝒮_I/N_err, 
                where i=1,2,......N_err;  

 
 
 (5) Variances: VA_L = 𝒮𝒮_L/N_err - ℳ_L,
 VA_M = 𝒮𝒮_M/N_err - ℳ_M,
 VA_R = 𝒮𝒮_R/N_err - ℳ_R,
 VA_S = 𝒮𝒮_S/N_err - ℳ_S,  and
 VA_I = 𝒮𝒮_I/N_err - ℳ_I;   

 
 
 (6) Co-variances:  CV_ML = 𝒮𝒮𝒮_ML/N_err - ℳ_M ×ℳ_L,
 CV_MR = 𝒮𝒮𝒮_MR/N_err - ℳ_M ×ℳ_R,
 CV_MS = 𝒮𝒮𝒮_MS/N_err - ℳ_M ×ℳ_S,
 CV_MI = 𝒮𝒮𝒮_MI/N_err - ℳ_M ×ℳ_I,
 CV_LR = 𝒮𝒮𝒮_LR/N_err - ℳ_L ×ℳ_R,
 CV_LS = 𝒮𝒮𝒮_LS/N_err - ℳ_L ×ℳ_S,
 CV_LI = 𝒮𝒮𝒮_LI/N_err - ℳ_L ×ℳ_I,
 CV_RS = 𝒮𝒮𝒮_RS/N_err - ℳ_R ×ℳ_S,
 CV_RI = 𝒮𝒮𝒮_RI/N_err - ℳ_R ×ℳ_I,
 CV_SI = 𝒮𝒮𝒮_SI/N_err - ℳ_S ×ℳ_I.
 


 §.§ Derivatives of β and log(L'_0) with respect to  the parameters M_s,  L, I_e, R_e,  and σ 


We start from relations (<ref>) and  (<ref>), however recast in the compact notation keeping the same units, and in which we have introduced  c_a = log(G/k_v), c_b = log(2π), and c_d = c_a + c_b. For the purposes of this analysis we keep constant the structural parameter k_v. Relations  (<ref>) and  (<ref>) become
	   

    log( L'_0)    =    
                [ (-4 R - 2 c_b) S -(I + M - L + R + c_d)×
           (2 S - M - I - R - c_d) ] / 
           [ (I  + M - L + R + c_d) - 2 S ]



    β   =    [1/(I + M - L + R + c_d)]×
          { [ 2 (4 R + 2 c_b) S + 2 (I + M - L + R + c_d) ×
                (2 S - M - I - R - c_d) ]/       
            [ (I + M - L + R + c_d) - 2S ] + (2 c_b + 4 R)  }
 
 
Secondly, we calculate the partial derivatives of log( L'_0) and β with respect to 
  log(M_s), log(R_e), log(σ), log(L),  and log(I_e), using the same compact notation. The derivatives are

  
    ∂log L'_0/∂ M   =      (2 c_b S + c_d^2 + 2 c_d (M + R - 2 S - L + I) + 
               M^2 + 2 M (R - 2 S - L + I) + R^2 - 2 R L + 2 R I + 
               4 S^2 + 2 S L - 4 S I + L^2 - 2 L I + I^2)/  
               ( c_d + M + R - 2 S - L + I)^2


  
    ∂log L'_0/∂ R   =    (2 c_b S + c_d^2 + 2 c_d (M + R - 4 S - L + I) + 
               M^2 + 2 M (R - 4 S - L + I) + R^2 - 4 R S - 2 R L +           
               2 R I + 12 S^2 + 6 S L - 8 S I + L^2 - 2 L I + I^2)/          
               (c_d + M + R - 2 S - L + I)^2

        

    ∂log L'_0/∂ S   =    -(2 (c_b + 2 R - L) (c_d + M + R - L + I))/  
              (c_d + M + R - 2 S - L + I)^2
	  


    ∂log L'_0/∂ L   =    (2 S (-c_b + c_d + M - R - 2 S + I))/         
                 (c_d + M + R - 2 S - L + I)^2



    ∂log L'_0/∂ I   =    (2 c_b S + c_d^2 + 2 c_d (M + R - 2 S - L + I) + 
              M^2 + 2 M (R - 2 S - L + I) + R^2 - 2 R L + 2 R I +      
              4 S^2 + 2 S L - 4 S I + L^2 - 2 L I + I^2)/              
              (c_d + M + R - 2 S - L + I)^2


and


    ∂β/∂ M   =     -(2 (c_b + 2 R - L))/(c_d + M + R - 2 S - L + I)^2



    ∂β/∂ R   =    -(2 (c_b - 2 c_d - 2 M + 4 S + L - 2 I))/   
               (c_d + M + R - 2 S - L + I)^2



    ∂β/∂ S   =    (4 (c_b + 2 R - L))/(c_d + M + R - 2 S - L + I)^2

	

    ∂β/∂ L   =    -(2 (-c_b + c_d + M - R - 2 S + I))/   
                (c_d + M + R - 2 S - L + I)^2



    ∂β/∂ I   =     -(2 (c_b + 2 R - L))/(c_d + M + R - 2 S - L + I)^2


Although it is not indicated in the notation, each partial derivative with respect to some variable  is calculated keeping constant all other variables. The total variations  of log L'_0 and β are expressed by 


    Δlog L'_0    =    (∂log L'_0/∂ M) Δ M + 
                 (∂log L'_0/∂ R) Δ R + 
                 (∂log L'_0/∂ S) Δ S +   
              (∂log L'_0/∂ L) Δ L + 
                 (∂log L'_0/∂ I)  Δ I



    Δβ   =    (∂β/∂ M) Δ M + 
                 (∂β/∂ R) Δ R + 
                 (∂β/∂ S) Δ S +   
              (∂β/∂ L) Δ L + 
                 (∂β/∂ I)  Δ I


These relations are derived for each object of the sample (i=1,2,.....N_max), and for each object i for all perturbations of its parameter set to simulate the uncertainty (j=1,2,.... N_err),

Particular care is paid to evaluate the variations Δ M, Δ L, Δ R, Δ S, and Δ to be inserted in relations (<ref> and (<ref>).
Since they all are written in logarithmic scale, their absolute value  cannot be larger than the uncertain θ=0.04. In reality it can be significantly smaller than this, depending on the method used to perturb the original values,  either uniform-random or gaussian-random distributions.



 §.§ Variance and co-variance of log L'_0 and β


If for the same galaxy we had many measurements of their characterizing parameters and an estimate of their errors, we could apply the technique of variances and covariances for correlated variables and  obtain an estimate of the propagation of errors on L'_0 and β, the target quantities of our analyses. To this aim we introduced the simulations of errors affecting the parameters (M, L, R, S, I, compact notation is used) and derived the associated variances and co-variances for each characterizing parameter of  each object in the sample. Thanks to it variance and co-variances of  log L'_0 and β of each galaxy and its simulations can be calculated. They are


    (VA_log L'_0)^2    =   ∑_J (∂log L'_0/∂ J)^2 VA_J +  
             
              +2 ∑_J,K(∂log L'_0/∂ J)(∂log L'_0/∂ K) CV_J,K



    (VA_β)^2    =   ∑_J (∂β/∂ J)^2 VA_J +  
             
              +2 ∑_J,K(∂β/∂ J)(∂β/∂ K) CV_J,K


where J=M,L,R,S,I and K=M,L,R,S,I, J ≠ K, all terms with J,K and K,J are equal and must counted only once. The terms  
VA_M , VA_L, VA_R, VA_S, VA_I are the variances,       and CV_ML,  CV_MR, CV_MS,  CV_MI, CV_LR, CV_LS, CV_LI, CV_RS, CV_RI, CV_SI  the co-variances. 
Finally, 

    (VA_log L'_0) = √((VA_log L'_0)^2)     and    
              (VA_β) = √((VA_β)^2)


Lastly, the percentage change  ΔΠ of β and log L'_o are given by


    ΔΠ_log L'_0   =   ((abs(Δlog L'_0 + log L'_0)-abs(log L'0))/(abs(log_L'0)100  
    ΔΠ_β   =   ((abs(Δβ + β)-abs(β))/(abs(β))100

	    


 §.§ Summary of the whole procedure


Given a galaxy with characterizing parameters,  M  = log(M_s), R  = log(R_e), S  = log(σ),  L  = log(L), I  = log(I_e), we  calculate 20 sub-cases whose characterizing parameters are artificially changed by random amounts always smaller than the estimated maximum uncertainty of 20% or 0.08 in logarithmic scale, that is ± 0.04 the unperturbed values. The generic parameters are here indicated by X where X stands for  M, R, S, L,  and I. For these moked objects we calculate the quantities: mean values 
    ℳ_X =  𝒮𝒮_X/N_err,
    variances
    VA_X = 𝒮𝒮_X/N_err - ℳ_X,    
    and covariances
    CV_X X' = 𝒮𝒮𝒮_X X'/N_err - ℳ_X ×ℳ_X',

where X' has the same meaning of X', the combinations  X ≠ X' are retained, and the combinations X X' = X' X are counted only once. Then we derived for the original and moked objects the formal solutions of the system of equations giving  log L'_0 and β and their partial derivatives, eqns. (<ref>) through eqn.(<ref>) and eqns.(<ref>) through (<ref>). With the aid of these, we derive the uncertainties affecting log L'_0 and β  and the propagation of uncertainties (errors in the estimates) for correlated variables. Finally, we derive the variances and co-variances of the solutions log L'_0 and β, eqns (<ref>). Lastly,  we calculate the percent variations of log L'_0 and β. The results were presented and commented in Sect. <ref>.
 
