© The Authors 2025

1. Introduction

Elliptical galaxies exhibit a well-known empirical relationship between their effective radius (R_e), central velocity dispersion (σ_e), and effective surface brightness (I_e), referred to as the fundamental plane (FP) (Djorgovski & Davis 1987; Dressler et al. 1987). The FP can be understood as a consequence of the virial equilibrium with homology and a constant mass-to-light ratio (M/L) among galaxies (Binney & Tremaine 2011). This relationship can be expressed as R_e ∝ σ_e^bI_e^c. According to the virial theorem, the theoretical values of these coefficients are b = 2 and c = −1. However, observational studies at low redshifts have empirically determined the values to be b ∼ 1.2 and c ∼ −0.8 (Jørgensen et al. 1996). Thus, observations reveal a tilt in the FP (Scodeggio et al. 1998; Trujillo et al. 2004; Hyde & Bernardi 2009; Cappellari et al. 2013) compared to the theoretical prediction for the virial equilibrium.

The evolutionary processes in elliptical galaxies, including changes in their structure, stellar populations, and mass-to-light ratios, cause FP slope parameters to evolve with redshift, leading to a rotation of the plane over cosmic time. In practice, values for the FP coefficients have been found to lie in the range b ≃ [0.4, 1.5] and c ≃ [ − 0.9, −0.4], depending on the chosen wavelength band and fitting method, as well as on redshift (van Dokkum & Franx 1996; Hyde & Bernardi 2009; Saglia et al. 2010; Jørgensen & Chiboucas 2013; Hou & Wang 2015; Stockmann et al. 2021). However, the evolution of the FP through observations has only been studied up to relatively low redshifts, z ∼ 2 (van de Sande et al. 2014; Stockmann et al. 2021). Spectroscopic observations needed to obtain high signal-to-noise spectra and measure stellar velocity dispersion at redshifts beyond 2 have proven challenging. So far, only the brightest galaxies have been studied at these early epochs (van de Sande et al. 2013; Belli et al. 2017; Stockmann et al. 2019), which might have resulted in a selection bias.

Most FP studies to date have focused on the luminosity-based FP (L-FP), which is highly sensitive to variations in stellar populations, dust attenuation, and observational wavelength. This becomes especially problematic at high redshifts, where rest-frame optical light shifts into the ultraviolet, tracing younger and more massive stellar populations and thereby failing to reliably trace the underlying mass distribution. The availability of stellar mass estimates from broadband photometry allows one to study the stellar mass FP (SM-FP), whereby surface brightness is replaced by the stellar surface mass density, Σ_e = M_*/2πR_e² (Hyde & Bernardi 2009; Bezanson et al. 2013; Hou & Wang 2016; Bernardi et al. 2020). The SM-FP is closer to the virial plane and enables a more physically meaningful interpretation of FP evolution.

Despite extensive observational efforts, the exact nature, origin, and evolution of the FP over cosmic time remain topics of active research. In recent years, various advanced cosmological hydrodynamical simulations including Horizon simulations (Dubois et al. 2014), EAGLE simulations (Schaye et al. 2015; Crain et al. 2015), the Next Generation Illustris (IllustrisTNG) simulations (Pillepich et al. 2018; Springel et al. 2018; Naiman et al. 2018; Marinacci et al. 2018; Nelson et al. 2018), the Magneticum Pathfinder simulation suite (Hirschmann et al. 2014; Teklu et al. 2015), and the Horizon Run 5 simulation (Lee et al. 2021; Park et al. 2022), incorporating complex physical processes such as star formation, feedback mechanisms, and interactions with dark matter provide a complementary approach to study the formation and evolution of galaxies (Snyder et al. 2015; Dubois et al. 2016; Correa et al. 2017; Park et al. 2022).

Recently, Lu et al. (2020) examined the L-FP (slopes, scatter, and zero-point) in the IllustrisTNG-100 simulation, finding a plane with mild evolution in its coefficients and a mass–plane relation in place, with a small scatter (∼0.08 dex) since z = 2. Their trends shows good agreement with observations at lower redshifts (z < 1), while at higher redshifts (z > 1) discrepancies arise, potentially due to the larger sizes (Genel et al. 2018; Rodriguez-Gomez et al. 2019) and lower σ_e (Wang et al. 2020) of IllustrisTNG early-type galaxies (ETGs) compared with observations.

Rosito et al. (2021) found that Horizon-AGN ETGs identified at z = 0 follow the observed fundamental scaling relations (mass–plane, FP, and mass–size relation). However, in simulations without AGN feedback, galaxies deviate from the expected FP, with some showing excess stellar surface density, leading them to the conclusion that active galactic nucleus (AGN) feedback plays a critical role in regulating galaxy size, rotation, and star formation activity, impacting the overall alignment with observed relations. This study also highlights how AGN feedback influences the mass-to-light ratio (M/L) and the tilt of the L-FP, noting a redshift-dependent evolution of the FP based on stellar age and morphology. However, most simulations have only extended FP studies up to a redshift of 2 ∼ 3 due to their small volume and limited particle resolution. Also, the definition of the ETGs is different among simulations and observations, making quantitative comparison of results difficult. Some recent works have provided important complementary perspectives on the FP across cosmic time. The Magneticum Pathfinder simulations (Dolag et al. 2025) analysed the SM-FP from z ∼ 4 to z = 0, reporting stable slopes, a steadily evolving zero point, and a constant intrinsic scatter in log R of ∼0.06 dex over the full redshift range, despite significant evolution in galaxy sizes and surface mass densities. In contrast, D’Onofrio & Chiosi (2024) investigated the luminosity- FP for both simulations (Illustris, IllustrisTNG) and observations, introducing a time-dependent scaling to connect FP tilt to the degree of virialization. Their framework predicts that bright massive ETGs remain close to the virial plane throughout cosmic history, while low-mass and low-luminous systems deviate more strongly, especially at high redshifts (z ∼ 4). The HR5 analysis presented here builds on these results by providing an independent measurement of the SM-FP coefficients and their redshift evolution.

This paper aims to study the redshift evolution of SM-FP of a class of galaxies selected based on their spheroidal morphology and/or internal kinematics by utilizing the recent Horizon Run 5 cosmological simulation (Lee et al. 2021) – HR5 hereafter – which uniquely combines a large volume and high resolution to probe galactic objects up to z ∼ 6. By adopting stellar mass surface density in place of luminosity-based surface brightness (also presented in Appendix A), we construct a SM-FP that is less sensitive to stellar population variations and more directly reflects the underlying dynamical state of galaxies. In Sect. 2, we describe the HR5 simulation and the selection criteria used to identify spheroidal and slow-rotating galaxies based on their morphology and kinematics. In Sect. 3, we outline the FP fitting methodology, including the use of robust least trimmed squares (LTS) plane fitting. In Sect. 4, we present our main results on the redshift evolution of FP coefficients for both a massive sub-sample and a broader sample of spheroids with a smaller mass cut, as well as the impact of a V_rot/σ_e-based kinematic selection criterion on the FP coefficients. In Sect. 5, we present a discussion emphasizing the importance of studying the SM-FP instead of the L-FP. We examine the mass dependence of the FP parameters, analyse how various stellar properties – such as the age, rotational velocity to stellar velocity dispersion ratio (V_rot/σ_e), and sSFR – vary along the FP, and compare our results with previous studies, including both observational data and simulations. We summarize our main conclusions in Sect. 6, highlighting the early emergence and persistence of the FP.

2. Data and methods

2.1. Simulation

The HR5 simulation (Lee et al. 2021) is a state-of-the-art cosmological hydrodynamical simulation that covers a volume of about (1048.6 cMpc)³, while achieving a resolution of ∼1 pkpc in the zoom-in region. The cosmological parameters adopted are based on the Planck observation (Planck Collaboration XIII 2016), with Ω_m = 0.3, Ω_Λ = 0.7, Ω_b = 0.047, σ₈ = 0.816, and h = 0.684. The initial conditions for the simulation are generated at z = 200 with the MUSIC package (Hahn & Abel 2011) using second-order Lagrangian perturbation theory (2LPT; Scoccimarro 1998; Jenkins 2010).

The HR5 simulation was run using the adaptive mesh refinement code RAMSES (Teyssier 2002). The RAMSES code was modified to implement OpenMP parallelism in addition to the Message Passing Interface already compatible with RAMSES. The Poisson equations and Euler equations were solved based on the particle mesh algorithm (Guillet & Teyssier 2011) and the Hartex-Lax-van Leer contact wave Riemann solver (Toro et al. 1994), respectively. All particles of dark matter, stars, and massive black holes were assigned into grids using the cloud-in-cell mass assignment scheme to construct mass-density grids.

The initial density field in a periodic box of 1048.6 cMpc size was generated on a coarse 256³ grid in which each cell size was 4.096 cMpc. A high-resolution cuboid region of volume 1048.6 × 119.0 × 127.2 cMpc³ was set at the centre of the simulation box. This high-resolution region has 8192 grids (dx = 0.128cMpc) along the long side in the IC. The simulation sub-grid physics parameters were set by specifically matching the cosmic star formation history (CSFH) with various observations. The detailed information about models of gas cooling and heating, AGN and supernova (SN) feedback, star formation, and chemistry are given in Lee et al. (2021). The galaxies in HR5 were identified from the entire snapshots of HR5 using the physically self-bound (PSB)-based galaxy finder, PGalF (Kim et al. 2023). PGalF uses the stellar density map to identify star particle groups (substructures or galaxies).

The strength of the HR5 simulation lies in the fact that it is able to capture physical processes on scales ranging from gigaparsec down to kiloparsec scales with a dynamic range an order of magnitude higher than any previous simulations of this kind. Figure 1 demonstrates the importance of using a simulation box larger than 1 cGpc. Plotted are the velocity two-point correlation functions (VCF) calculated from the Fourier transformation of the power spectrum P(k)/k², integrated over the ranges covered in the initial conditions of recent simulations such as Horizon Run 4 (green; 3150 h⁻¹Mpc; Kim et al. 2015), HR5 (red; 717.2 h⁻¹Mpc), Illustris TNG300 (brown; 205 h⁻¹Mpc), TNG100 (violet; 75 h⁻¹Mpc), and TNG50 (pink; 35 h⁻¹Mpc; Nelson et al. 2019). It is related with the pairwise peculiar velocity v₁₂(r) by ξ_v(r) = ξ_v(0)−⟨v₁₂²⟩/2. The black line in the upper panel is the correct VCF when the full spectral range is covered, and the VCFs relative to the correct one endowed in the initial conditions of various simulations are shown in the bottom panel. It can be seen that the velocity dispersion, or the VCF at zero separation, is significantly underestimated when the simulation box size is smaller than 1 cGpc. This may be one of the reasons why the galaxy velocity dispersion is on average lower in TNG50 compared to that in TNG100 at a fixed galaxy mass even though TNG50 has a higher force resolution than TNG100. The same is true between TNG100 and TNG300. At the galaxy (∼1 cMpc) and cluster (∼10 cMpc) scales, the VCF is even more underestimated, and it is hard to anticipate a correct calculation of the internal kinematic evolution of galaxies and clusters in the simulations with box size much smaller than 1 cGpc as mergers and accretion are not accurately simulated.

Fig. 1. Comparison of VCF for the various cosmological hydrodynamical simulations highlighting the importance of larger simulation box size. It can be seen that the VCF is significantly underestimated at zero separation (velocity dispersion), galaxy scale (r ≈ 1 h−1 Mpc) and cluster scale (r ≈ 10 h−1 Mpc) when the simulation size is less than 1cGpc

Another critical benefit of using HR5 is that it allows us to examine the FP at high redshifts. Due to its large volume, a significant number of galaxies have already formed by z = 6, and about 15% of them have a spheroidal morphology (Park et al. 2022). This allows us to study the FP of spheroids at least up to z = 6. However, it should also be pointed out that the force resolution of HR5 is about 1 pkpc and the half mass radius of small galaxies could have been systematically affected. Therefore, we aim to study the emergence of the FP and relative evolution of the FP parameters rather than their absolute values.

2.2. Sample selection

In this section we outline our criteria for selecting the sub-sample of galaxies from the HR5 simulation for our FP study. Initially, we applied a threshold on the galaxy stellar mass, considering only central galaxies with stellar mass M_* > 2 × 10⁹ M_⊙ or 1,000 stellar particles, to securely estimate the internal structure of galaxies and to reduce the potential impact of the limited resolution of the simulation.

Following the methodology of Park et al. (2022), we further classified the galaxy morphology based on the asymmetry parameter and the radial profile of their stellar mass density distribution. We utilized a non-parametric morphological descriptor, the asymmetry parameter (A), to distinguish between symmetric and asymmetric galaxies. Park et al. (2022) observed a dichotomy in A at approximately A = 0.4 to 0.5 for galaxies with a relatively high stellar mass. Accordingly, galaxies with A ≥ 0.4 are deemed to be asymmetric, while those with (A < 0.4) are considered symmetric.

Subsequently, the target sample was refined by fitting the radial profile of each galaxy’s stellar mass density projected on the galactic plane with a Sérsic profile, allowing us to derive the best-fit Sérsic index value, n_Sérsic. The regular galaxies with A < 0.4, 1.5 < n_Sérsic < 100 are classified as spheroids, and those with 0 < n_Sérsic < 1.5 are called the disc galaxies. The value of the n_Sérsic = 1.5 chosen to classify galaxies into two broad morphology classes is based on the distribution of the Sérsic index of SDSS galaxies (Park & Choi 2005). The radial range over which the stellar mass-density profile was computed was constrained by physical and numerical considerations. R₉₀, containing 90% of total galaxy stellar mass, is the outer limit of the radial range where the mass-density profile was calculated and used to estimate the Sérsic index. The 10% of stellar particles at the outermost region were dropped to minimize the mass contribution from possible stray galaxies and the tidal disruption effects of neighboring systems. On the other hand, the inner limit was set to r = 0.8 pkpc, taking into account the force resolution of HR5 and to avoid fitting galaxy cores. The above galaxy morphology classification criteria reproduce the fraction of discs, spheroids, and irregular galaxies as a function of stellar mass and redshift, which is in close agreement with the observational results from JWST (Lee et al. 2024).

The lowest redshift considered in our study is z = 0.625, which corresponds to the final snapshot available from the HR5 simulation that evolves from z = 200 down to z = 0.625. In Fig. 2, we present key properties such as the half stellar mass radius (R_hsm), SDSS (g − r) colour, SDSS r-band age (age_r), and Sérsic index (n_Sérsic) versus stellar mass (M_*) for our sample selected at z = 0.625. The blue points represent the spheroids sample with M * > 2 × 10⁹ M_⊙, while the orange points show the distribution of the spheroid sub-sample with M * > 10¹¹ M_⊙. Tables 1 and 2 summarize the number of spheroids in each sub-sample at each redshift analysed in our analysis. In addition to the morphology-based samples, we also constructed sub-samples using an additional criterion of V_rot/σ_e < 0.55 to identify slow rotators, thereby examining the impact of internal kinematics on the FP coefficients. This threshold is consistent with previous studies (Chisari et al. 2015; Rosito et al. 2021) and, in our case, also corresponds to the median V_rot/σ_e value measured for HR5 galaxies as a function of stellar mass. The rotational velocity, V_rot, is the mass weighted tangential velocity of stellar particles inside the 3D stellar half mass radius in cylindrical co-ordinates. The velocity dispersion, σ_e, is defined as ${\sigma }_e=\sqrt{({\sigma }_R^2+{\sigma }_{\mathrm{rot}}^2+{\sigma }_z^2)/3}$, where σ_R, σ_rot, and σ_z represent the velocity dispersion in the radial, tangential, and vertical directions, respectively.

Fig. 2. Distribution of 3D half stellar mass radius (log Rhsm), SDSS (g − r) colour, stellar age in r band (log ager), and Sérsic index (nSérsic) versus stellar mass (log M*) for selected spheroidal galaxies (blue points), their massive sub-sample (${\mathrm{M}}_{\ast }>{10}^{11}{\mathrm{M}}_{\odot }$; orange points), and the full sample of galaxies (with ${\mathrm{M}}_{\ast }>2\times {10}^9{\mathrm{M}}_{\odot }$; black contours) from the HR5 simulations at z = 0.625, respectively. The dashed black line in the bottom right-hand panel indicates nSérsic = 1.5.

3. FP fitting

The observed FP relation is given by

$$\begin{array}{c}\hfill log{R}_{\mathrm{e}}=a+blog{\sigma }_e+clog⟨I_e⟩,\end{array}$$(1)

where R_e is the effective radius, σ_e is the 1D line of sight velocity dispersion measured from the galaxy spectra, and ⟨I_e⟩ is the mean surface brightness within R_e defined as I_e = L_e/2πR_e².

To resolve the contribution of the stellar population, the SM-FP can be measured by replacing the luminosity term in the above L-FP with stellar mass, such that

$$\begin{array}{c}\hfill log{R}_{\mathrm{hsm}}=a+blog{\sigma }_e+clog{\mathrm{\Sigma }}_e,\end{array}$$(2)

where Σ_e = M_*/(2πR_hsm²) is the mean stellar mass density, with M_* being the total stellar mass within the half-stellar-mass radius, R_hsm (measured in solar masses), σ_e the mass-weighted root-mean-square velocity of all stellar particles within R_hsm calculated along a given line of sight (z axis, measured in kilometres per second), and R_hsm is the 3D half stellar mass radius (measured in kiloparsecs). We then used the LTS-Plane fit (Rousseeuw & Van Driessen 2006; Cappellari et al. 2013) method to perform mass FP plane fitting on our galaxies at various redshifts. This method uses the least-squares fitting algorithm and also trims away outliers that contribute most to the residuals. For each analysed redshift, we estimated the parameters defined in the Eq. (2) for spheroids in our sample.

To estimate the uncertainties of the FP coefficients, we adopted two independent methods. In the first approach, we generated 40 bootstrap resamples, each with the same size as the original dataset, and adopted the standard deviation of the fitted FP coefficients across these resamples as the parameter uncertainties. In the second approach, we divided the sample at each redshift into n = 10 sub-samples of equal size, fitted the FP independently to each, and calculated the standard deviation of the resulting coefficients. This procedure was repeated ten times, and the average of the standard deviations was taken as the final uncertainty estimate. A comparison of the two approaches shows that the bootstrap errors are systematically smaller than those from sub-sampling, indicating that bootstrapping tends to underestimate the true uncertainties in the FP coefficients. Nevertheless, we report results from both methods and specify which one is adopted in each case. We present our results in the next section.

4. Results

We now present the best fitting FP parameters, b and c, its zero point, a, as well as the scatter, Δ, and relative scatter, Δ/Δ_r, which can be calculated as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Delta }& =\sqrt{\frac{\sum _i^N{(a+blog{\sigma }_{\mathrm{e},\mathrm{i}}+clog{I}_{\mathrm{e},\mathrm{i}}-log{R}_{\mathrm{hsm},i})}^2}{N}},\hfill \\ \hfill {\mathrm{\Delta }}_r& =log{R}_{\max }-log{R}_{\min },\hfill \end{array}\end{array}$$(3)

where N is the total number of selected galaxies at each redshift, and R_max and R_min mark the range of R_hsm that contains 95% of the total sample of galaxies computed using 2.5 and 97.5 percentile values, respectively. We chose Δ_r to estimate the range in log R_hsm instead of more commonly used σ_R, which is the RMS value. Our choice of Δ_r was motivated by the following reasons: (i) it directly quantifies the extent of galaxy size distribution in the sample, offering a clearer picture of the structural diversity, and (ii) it is less sensitive to outliers or to the strongly skewed distribution of galaxy radii, whereas σ_R undesirably decreases when small-size galaxies dominate the sample.

4.1. Stellar mass fundamental plane: Massive spheroids with M_* > 1 × 10¹¹ M_⊙

We present in the left panel of Fig. 3 the redshift evolution of the FP of spheroidal galaxies with stellar masses of M_* > 10¹¹ M_⊙ across various redshifts. The simulated FPs are well defined, exhibiting a relatively low scatter at all redshifts up to z ≃ 2.5. We observe a shift in the FPs towards the right as the redshift decreases, indicating the increasing size evolution of galaxies over cosmic time. Numerous studies have shown that at a fixed stellar mass, spheroid-like galaxies were significantly more compact at higher redshifts (z ∼ 1 − 2) than their likely present-day descendants (Daddi et al. 2005; Longhetti et al. 2007; Trujillo et al. 2007; Buitrago et al. 2008; Stockmann et al. 2019, 2021). The drivers of spheroid galaxy size evolution in our simulations will be explored in a separate paper.

Fig. 3. Redshift evolution of SM-FP for the spheroidal galaxies in HR5. Note that galaxies are located along the diagonal line as the best-fit FP coefficients are used at each redshift. The relation moves to the right as galaxies grow in size with time. Left: For the massive sub-sample with M* > 1 × 1011 M⊙Right: For the sample with M* > 2 × 109 M⊙. To mitigate the overlap of points, black contours are plotted for the initial five redshifts, representing the central 90% of galaxies contained within.

Table 1 summarizes the FP slope parameters, b and c, the zero point, a, the absolute scatter, Δ, and the relative scatter, Δ/Δ_r, for the massive galaxies. The error bars in this case were obtained using bootstrap resampling method as is described in Sect. 3 because of small sample sizes.

To examine the evolution of the FP slope parameters, and to facilitate comparisons with other simulations and observational studies, we present in Fig. 4, the redshift dependence of coefficients, b (upper panel) and c (middle panel). We find that in the case of our massive sub-sample, the overall magnitude of the coefficient, c, represented by the solid red curve is smaller than the Horizon-AGN (Rosito et al. 2021) result represented by the dashed green curve, at all redshifts. However, both simulations predict a weak (possibly decreasing) dependence of |c| towards higher redshifts. It should be pointed out that the sample selection criteria for the FP study are different in different simulations, and hence we do not expect to see an accurate quantitative agreement among the simulated results. The Horizon-AGN study (Rosito et al. 2021) distinguishes slow rotators from fast rotators, using V_rot/σ_e to classify them into early or late types, respectively, and follows Chisari et al. (2015) in adopting V_rot/σ_e < 0.55 to select a slowly rotating galaxy sample. In contrast, our methodology selects spheroidal galaxies based directly on their morphology using structural parameters, including the asymmetry parameter and the Sérsic index. These differing selection criteria may introduce systematic differences in the galaxy populations analysed, thereby influencing the quantitative comparison of FP parameter evolution between simulations.

Fig. 4. Redshift evolution of the SM-FP parameters: slope coefficients, b (top) and c (middle), and zero-point variation, Δa = az − azi (bottom). The results of the HR5 samples are indicated by a red curve for the massive sub-sample and a blue curve for the full sample. The errors on a, b, c, Δ, and Δ/Δr were estimated using bootstrap resampling for the massive sub-sample, while for the broader sample that includes lower-mass galaxies, they were determined with the sub-sampling method. Observational constraints from Hou & Wang (2016) and Bernardi et al. (2020) are shown as magenta and yellow points, respectively.

On the other hand, the magnitude and redshift evolution of coefficient b agrees quite well with the Horizon-AGN predicted evolution. An interesting finding is that both simulations predict a weak dependence of b on redshift in the sense that b monotonically decreases toward higher redshifts. We show below that this trend of the b coefficient is still true for less massive galaxies but with a different magnitude of b. This indicates that the coefficient b reflects galaxy evolution more sensitively than c. Similarly, the FP zero point offset defined as Δa = a_z − a_zi, represented by the solid red curve in the bottom panel of Fig. 4, tracks the FP zero-point shift, isolating any structural or dynamical changes at fixed σ_e and Σ_e. We find negligible evolution of Δa, suggesting early structural maturity for the massive galaxies. Another interesting property of the FP is its tightness (i.e., the amount of scatter of galaxies around the plane). It is related to questions such as how early an FP forms, and how its tightness evolves with cosmic time sheds light on galaxy-evolution processes. We calculated the absolute scatter, Δ, to quantify the strength of the FP and the relative scatter, Δ/Δ_r, to quantify the existence of the FP. We note from Table 1 that both Δ and Δ_r remain low at all analysed redshifts, indicating the existence of a well-defined FP at z ≃ 2.5 (see the left panel of Fig. 3). This in turn suggests that massive galaxies were already structurally mature and dynamically settled at this epoch, so that the FP was established very quickly in cosmic history.

The FP parameters measured for observed galaxies are too scarce for a direct comparison with simulations, especially at redshifts higher than z ∼ 1. Additionally, the comparison is also impacted by the fact that the observational results have used the ETG samples selected using different criteria. For instance, in a FP study based on SDSS DR7 galaxies at z < 0.2 (Hou & Wang 2015), the ETGs were selected based on multiple criteria, i.e. (a) the concentration index r₅₀/r₉₀ > 2.5, (b) the ratio of the likelihood of the de Vaucouleurs model to the exponential model L_deV/L_exp ≥ 1.03, and (c) a spectral classification index eClass of less than −0.1.

4.2. Stellar mass fundamental plane: Spheroids with M_* > 2 × 10⁹ M_⊙

Our full sample contains the spheroidal HR5 galaxies with M_* > 2 × 10⁹ M_⊙. The right panel of Fig. 3 shows that the FP remains well defined across all redshifts, from z = 0.625 up to z = 6, with a scatter that is moderate but larger than that of massive ones. At z ∼ 6, most galaxies have M_* < 1 × 10¹⁰ M_⊙, and only about 15% of them have a spheroidal morphology. These spheroidal galaxies tend to undergo morphological transformation back to disc types as the large-scale tidal torques set up in the primordial density fluctuations continue to inject angular momentum via subsequent mergers and accretion (Park et al. 2022). Even though the spheroidal morphology of these galaxies is incidental and temporary, we find that the spheroidal galaxies at high redshifts still lie on a well-defined FP. Therefore, the right panel highlights the robustness of the FP across a wide mass range. It shows how galaxy formation and evolution processes leave a clear imprint on the FP even at the earliest cosmic epochs.

Although the FP is already established by z ∼ 6, its normalization still shifts as galaxies grow. As redshift decreases, the FP shifts towards larger log R, indicating an overall increase in galaxy size over cosmic time, consistent with observations that galaxies tend to grow larger at later epochs. In Table 2, we summarize the FP coefficients, a, b, and c, the absolute scatter, Δ, and the relative scatter, Δ/Δ_r, across a range of redshifts, extending up to z ∼ 6. The uncertainties on each coefficient were estimated in the manner described in Sect. 3.

In Fig. 4, the blue curve represents the redshift evolution of the FP coefficients b and c, as well as the zero-point offset, Δa, for this broader galaxy sample. The magnitude of c in this case agrees (within uncertainties) with the value measured for massive galaxies (M_* > 10¹¹ M_⊙) up to z ∼ 2, though |c| remains larger than the Horizon-AGN prediction. Also, there is an indication of a decrease towards high redshifts beyond z > 1.2. Similarly to the massive ones, b decreases gradually as redshift increases, but its value is larger than in the HR5 massive sub-sample or the Horizon-AGN slow rotator track. Meanwhile, Δa, the FP zero-point offset, increases with redshift up to z ≈ 2, then decreases at higher z. This trend likely results from structural evolution at fixed σ_e and Σ_e, as well as from changes in the dynamical-to-stellar mass ratio (e.g. an evolving dark-matter fraction).

To further refine our sample, we examined the kinematics of the morphologically selected spheroids (A < 0.4, n_s > 1.5) and found that about 50% of them are highly dispersion-dominated (V_rot/σ_e < 0.55). We applied this additional kinematic selection to isolate dispersion-dominated spheroids. This criterion yields a subset of spheroids with a tighter FP, with the relative scatter decreasing by about 10–15% across the redshift range. The FP coefficients shift only slightly, bringing them closer to those of the massive galaxy sub-sample, as is indicated by black curves in each panel of Fig. 4. The overall redshift evolution of the FP coefficients remains qualitatively similar after applying this additional selection criterion.

5. Discussion

In the previous section, we reported the existence of the FP for spheroidal galaxies in the HR5 simulation at redshifts beyond 3 and as high as z ∼ 6. To our knowledge, no previous galaxy simulations or observations have explored the FP at such early cosmic epochs. Both the absolute scatter (Δ) and the relative scatter (Δ/Δ_r) remain low even at high redshifts, as is summarized in Table 2. The finding that the mass–FP is already well established by z ∼ 6 is striking, since galaxies at these epochs are still actively forming and undergoing frequent mergers. However, it should be noted that most galaxies have M_* < 10¹⁰ M_⊙ at high redshifts of z ∼ 6, and there are essentially no AGN activities or rapid star formation quenching active in HR5 simulation at these epochs. The existence of FP therefore reflects the internal dynamical stability of the stellar components of galaxies. A natural explanation lies in the processes that rapidly drive spheroidal systems towards equilibrium. In particular, violent relaxation during major mergers efficiently erases the memory of progenitor orbits and establishes virialized systems whose global properties are largely independent of the detailed merger history. At high redshift, dissipative (gas-rich) mergers further promote compact remnants, thereby tightening the scaling relations.

5.1. Reasons for studying the stellar mass FP instead of the luminosity FP

The conventional luminosity FP has long been used as a tool to study ETG evolution. However, the L-FP is sensitive to variations in stellar populations, dust attenuation, and observational wavelength (Scodeggio et al. 1998; La Barbera et al. 2008; Hyde & Bernardi 2009). This becomes especially problematic at high redshifts, where the spectral energy distribution of galaxies evolves rapidly, and even the rest-frame luminosity at a particular band does not reliably trace the underlying stellar mass distribution after the cosmological band-shift correction is made. To mitigate these biases, we instead focus on the SM-FP, replacing surface brightness, I_e, with the stellar mass surface density, Σ_e. This allows us to isolate structural and dynamical properties independently of the spectral energy distribution and contamination from dust. As is shown in Appendix A, the r-band luminosity becomes an increasingly poor tracer of stellar mass at z ≥ 1, whereas the SM-FP remains robust out to z ∼ 6. This mass-based approach provides a more direct probe of the underlying dynamically relaxed baryonic mass component and yields a more physically interpretable FP.

5.2. Mass dependence of FP parameters

Our analysis reveals a strong mass dependence in the FP slope, b, and in its zero-point offset, Δa, as is shown in Fig. 4 and Tables 1 and 2. First, the coefficient b is about 60% larger for the whole sample compared to that of the massive subset. Nevertheless, its redshift evolution does not seem to depend on mass, as both samples show a decreasing trend with redshift for b with similar slopes. On the other hand, the zero point offset, Δa, shows a strong dependence on mass in both its magnitude and evolution rate, as is shown in the bottom panel of Fig. 4. Interestingly, we find much less dependence of c on stellar mass over the common redshift range explored for the two mass samples.

5.3. Nature of the FP scatter

In addition to slope parameters, we find that the FP scatter in HR5 remains low (Δ ∼ 0.05 − 0.07 dex) for the whole sample and as well as (∼0.03 dex) for the massive sample, even at high redshifts, suggesting an early emergence of a well-defined mass–plane relation. This tightness is consistent with the small scatter (∼0.06 dex) found in the Horizon-AGN FP study up to z ∼ 3, although HR5 uniquely extends the FP analysis up to z ∼ 6, where such tight scaling relations have not previously been reported.

To study the nature of FP scatter, we examined the distribution of the stellar age, V_rot/σ_e ratio, and sSFR of galaxies along the FP (see Fig. 5). These properties are often used to select galaxy samples for an FP analysis. For instance, Lu et al. (2020) used the sSFR to select the quenched galaxies for the FP study in the Illustris-TNG simulation, while Rosito et al. (2021) used the V_rot/σ_e ratio to isolate slow rotators in the Horizon-AGN simulation. In both cases, the selected galaxies are broadly called the ‘early-type’ galaxies. In our study, we adopted the morphology of galaxy stellar mass distribution to select spheroidals for the FP analysis.

Fig. 5. Redshift evolution of SM-FP for the spheroidal galaxies in HR5. The symbols are coloured according to the mass-weighted average stellar age (upper panels), Vrot/σe (middle panels), and sSFR = SFR/M*[yr−1] (bottom panel). The 1:1 line is depicted in black in each panel.

At z ∼ 1, Fig. 5 shows a noticeable age gradient along the FP, with galaxies having older stellar populations systematically located towards larger effective radii and lying mostly above the 1 : 1 relation. Beyond z ≈ 4, this trend is reversed, with relatively older galaxies migrating below the 1 : 1 line and to smaller radii. This interesting transition happens starting from the cosmic noon period between z ∼ 4 and 1.5. A similar trend is visible in the sSFR shown in the bottom panels of Fig. 5, where the galaxies with a relatively low sSFR show a very similar behaviour. We have also found a similar trend along the FP for the stellar mass-to-light ratio (M_*/L). This is because older galaxies exhibit a systematically higher M_*/L ratio. These findings agree with the trends observed in the Horizon-AGN (Rosito et al. 2021) study, where AGN feedback was shown to drive older stellar ages and tighter FP relations.

The distribution of the age, sSFR, and M_*/L ratio along the FP indicates that galaxies are not cleanly separated by these properties. At z > 3, a subset of relatively passive or older galaxies defines a much tighter FP than the full sample. However, at z ≲ 3, a restriction in these parameters results in nontrivial evolution of the subset along the FP, and does not reduce its overall scatter. This implies that the age, sSFR, and M_*/L ratio are not the fundamental physical drivers of the FP.

On the other hand, the middle panels of Fig. 5 clearly show a systematic dependence of FP tightness on V_rot/σ_e. At all redshifts, slow rotators define a tighter FP compared to the whole sample and lie above the 1 : 1 relation. Below z ∼ 2, however, this distinction weakens as the slow rotator population loses some of its smaller members. This indicates that it is the V_rot/σ_e parameter that better selects the tighter FP galaxies compared to the star formation activity indicators. Constraining our sample of all spheroidal-type galaxies with M_* > 2 × 10⁹ M_⊙ to those with V_rot/σ_e < 0.55 results in a tighter FP with a smaller Δ (higher strength of FP), closer to that of the massive sample, as can be seen by the black curves in Fig. 4.

5.4. Comparison with other studies

Here, we compare our SM-FP results with predictions from other cosmological hydrodynamical simulations and some observational studies. While differences in sample selection and physical models prevent a one-to-one quantitative comparison, several key trends emerge.

For both HR5 and Horizon-AGN simulations, the magnitude of the FP coefficient, |c|, which quantifies the dependence on the stellar mass surface density, Σ_e, shows a mild decrease with redshift. The coefficient b, representing the dependence on the velocity dispersion, σ_e, exhibits a redshift evolution in both HR5 and Horizon AGN. For the HR5 full sample (including lower-mass galaxies), b decreases steadily from ∼2.1 at z ∼ 0.625 to below 1.0 at z ∼ 4 to 6, indicating a shallower relation between galaxy size and velocity dispersion at early times. This trend is consistent with the Horizon-AGN results, though the whole sample of HR5 galaxies with M_* > 2 × 10⁹ M_⊙ tends to show higher b values at fixed redshift. For the massive sub-sample with M_* > 10¹¹ M_⊙, the b values from HR5 are much closer to Horizon AGN results. Additional insights come from the Magneticum Pathfinder simulations (Dolag et al. 2025), which examined the SM-FP from z∼4 to z = 0. They find nearly constant slopes, a steadily increasing zero point towards higher redshift, and a constant intrinsic scatter of ∼0.06 dex in log R_e, despite strong size and surface density evolution. While differing observable definitions hinder a direct numerical comparison, both HR5 and Magneticum simulations show a tight SM-FP up to high redshift, with the evolution driven mainly by zero-point shifts rather than slope changes. The two observational studies (Hou & Wang 2016; Bernardi et al. 2020) of the SM-FP plotted in Fig. 4 are both at redshifts lower than the lowest HR5 snapshot (z = 0.625). A direct comparison with our results is therefore not possible, but they provide useful context for the FP at low redshift. Looking forwards, with the advent of JWST (NIRSpec/NIRCam; Jakobsen et al. 2022; Rieke et al. 2005) and large ground-based facilities such as Gemini/FLAMINGOS-2 (Eikenberry et al. 2006), it will become possible to measure high signal-to-noise velocity dispersions and structural parameters for galaxies at z > 2. These datasets will provide a critical test of our prediction that a tight SM-FP is already in place by z ∼ 6.

In Appendix A we also present the redshift evolution of the L-FP to compare with our SM-FP results. Figure A.1 illustrates the evolution of the L-FP coefficients b and c, as well as the zero-point offset, Δa, for both the massive galaxy sub-sample and the broader sample. We find that, consistent with previous simulation studies, the coefficient c exhibits negligible redshift dependence across all samples, while the coefficient b shows a systematic decrease with increasing redshift, indicating sensitivity to structural evolution. Observational results from van Dokkum & Franx (1996), Saglia et al. (2010), Hou & Wang (2015), and Stockmann et al. (2021) are also included for comparison. While these lie broadly within the simulation predictions at low redshift, the scarcity of high-redshift data limits a direct comparison beyond z ∼ 2.5. The zero-point offset, Δa, remains nearly flat for the massive galaxies, suggesting early structural maturity. In contrast, the lower-mass sample shows a pronounced increase in Δa up to z ∼ 2.5, followed by a decline, possibly reflecting transitions in the dominant physical processes such as mergers or evolving mass-to-light ratios. In Fig. 6, we compare the SM-FP (solid lines) with the corresponding L-FP (dashed lines). The differences are evident. For the coefficient c, the SM-FP evolves gradually and already approaches the theoretical virial expectation value of −1 by z ∼ 2, whereas the L-FP shows a stronger offset and slower convergence. A similar contrast is seen in the zero-point offset, Δa: the SM-FP evolves smoothly with a modest rise, while the L-FP exhibits a much larger variation, before declining towards zero at lower redshift. For the coefficient b, the SM-FP increases steadily and reaches higher values at z > 2, converging on the virial expectation of b ≈ 2. Although the SM-FP coefficients (b and c) vary more strongly with redshift than those of the L-FP, their evolution follows a clear, physically consistent trend towards the virial expectations (b ≈ 2, c ≈ –1). In contrast, the L-FP evolves more slowly and remains systematically offset from these values. This indicates that the SM-FP traces the structural and dynamical evolution of galaxies more coherently, whereas the L-FP retains the imprint of stellar population evolution.

Fig. 6. Comparison between the L-FP and the SM-FP as a function of redshift. The panels show the evolution of FP coefficients b (top) and c (middle), and the zero-point offset, Δa = az − azi (bottom). Results are shown for two stellar mass-selected samples: a massive sub-sample with M* > 1 × 1011 M⊙ (red) and the full sample with M* > 2 × 109 M⊙ (blue). Solid lines correspond to the SM-FP, while dashed lines indicate the L-FP.

6. Conclusions

We have investigated the FP of the spheroidal morphological type galaxies at redshifts up to z = 6, using the Horizon Run 5 simulation, which captures the initial density fluctuation modes up to ∼1cGpc. Thanks to its unprecedented large dynamic range from the gigaparsec scale to the scale of ∼1 proper kpc, HR5 allows us to identify a significant number of galaxies with M_* > 1 × 10⁹ M_⊙, even at redshift 6 and greater.

To examine the redshift evolution of the FP, in each HR5 snapshot we selected a subset of galaxies satisfying the following two conditions: (i) being the central galaxy, and (ii) having a spheroidal morphology (n_Sérsic > 1.5, where n_Sérsic is the Sérsic index of the stellar mass density profile at distances larger than 0.8 pkpc from the galaxy centre). It should be noted that we selected galaxies solely by their morphology and did not require them to be quiescent based on their colour or sSFR-properties that do not necessarily correlate with morphology or internal kinematics, particularly at early cosmic epochs. As the relations between the star formation rate, internal kinematics, and morphology evolve with time and the well-known tight relations are merely those near the present epoch, one should not use one of these properties as a proxy of the other at high redshifts Park et al. (2022).

In this study, we find that the FP relation was already well established at z ≈ 6, with a small scatter. Additionally, we find that the coefficient b of the log σ_e term exhibits a systematic decreasing trend with increasing redshift. The magnitude of the coefficient c of the Σ_e term shows a similar systematic decrease. These trends reflect the evolving structural and dynamical properties of spheroidal galaxies over cosmic time. The absolute values of the FP coefficients seem to depend on galaxy mass. Coefficient b is smaller for more massive galaxies with M_* > 10¹¹ M_⊙, compared to a broader sample of galaxies with M_* > 2 × 10⁹ M_⊙, while the two samples have similar values of coefficient c. Different simulations tend to give different absolute values of the coefficients, but interestingly predict more or less the same redshift dependence. That is, all simulations predict decreasing evolution of the magnitudes of b and |c| as redshift increases. On the other hand, a direct comparison with observations remains challenging due to differences in sample selection criteria and the limited availability of data at z > 2.

Our main findings are summarized as follows:

1. A well-defined SM-FP exists as early as z ∼ 6, with a low intrinsic scatter (Δ ≲ 0.07 dex), indicating that the dynamical and structural correlations of spheroidal morphological-type galaxies and slowly rotating galaxies were already established in the first billion years of cosmic history.

2. We find a mass-dependence of the FP coefficients. Coefficient b, representing the velocity dispersion term, is smaller for relatively massive galaxies with M_* > 10¹¹ M_⊙ compared to those for galaxies with M_* > 2 × 10⁹ M_⊙. The zero point, a, shows almost no evolution for the massive ones, while the smaller mass galaxies show non-monotonic evolution in a. In contrast, the coefficient c, tied to the stellar surface mass density, shows nearly no mass dependence.

3. The magnitudes of b (positive) and c (negative) show a roughly monotonic decreasing trend with redshift, while a for relatively low-mass galaxies increases up to z ≈ 2 and decreases above that.

4. Comparisons with other simulations, including Horizon-AGN and IllustrisTNG, show broad agreement in the redshift trends of FP coefficients, despite differences in galaxy selection criteria. HR5 uniquely enables the study of the FP at z > 2.5, where observational constraints are currently limited.

5. An analysis of FP residuals reveals systematic trends with the stellar age, V_rot/σ_e ratio, and sSFR (and also mass-to-light ratio), supporting the interpretation that stellar population properties influence the FP zero point and scatter, even in mass-based scaling relations. Among the parameters, V_rot/σ_e reduces the FP scatter at all redshifts.

The early emergence and persistence of the SM-FP can provide key constraints for galaxy formation models and motivate high-redshift spectroscopic observations with JWST and ground-based facilities such as Gemini to test these predictions.

Acknowledgments

PG and CP are supported by KIAS Individual Grant (PG088101, PG016903) at Korea Institute for Advanced Study (KIAS). CP is also supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT; RS-2024-00360385). J.L. is supported by the National Research Foundation (NRF) of Korea grant funded by the Korea government (MSIT, RS-2021-NR061998). S.E.H. is partly supported by the Korea Astronomy and Space Science Institute grant funded by the Korean government (MSIT; No.2024186901) and the grant funded by the Ministry of Science (No. 1711193044). This work benefited from the outstanding support provided by the KISTI National Supercomputing Center and its Nurion Supercomputer through the Grand Challenge Program (KSC-2018-CHA-0003, KSC-2019-CHA-0002). It was also supported by the Supercomputing Center/Korea Institute of Science and Technology Information, with supercomputing resources including technical support (KSC-2013-G2-003), and the simulation data were transferred through a high-speed network provided by KREONET/GLORIAD. This work also made use of computing resources at the Center for Advanced Computation at Korea Institute for Advanced Study. Software: AstroPy (Astropy Collaboration 2022), Matplotlib (Hunter 2007), Pandas (The pandas development team 2020; McKinney 2010), Seaborn (Waskom 2021).

References

Appendix A: Luminosity fundamental plane

We also investigate the conventional L-FP by adopting the following mathematical form:

$$\begin{array}{c}\hfill log\frac{R_{\mathrm{hsm}}}{[\mathrm{kpc}]}=a+blog\frac{{\sigma }_e}{[\mathrm{km}{\mathrm{s}}^{-1}]}+clog\frac{I_e}{[L_{\odot }{\mathrm{kpc}}^{-2}]},\end{array}$$(A.1)

where again σ_e is the mass-weighted root-mean-square velocity of all stellar particles within R_hsm; I_e is the rest-frame r-band mean surface brightness of all stellar particles within R_hsm, defined as follows:

$$\begin{array}{c}\hfill I_e=\frac{L}{\pi R_{\mathrm{hsm}}^2},\end{array}$$(A.2)

where L is the rest-frame r-band luminosity of the galaxy projected within R_hsm. The r-band luminosity is not an ideal tracer of stellar mass especially at high redshifts no matter how the early-type population of galaxies are defined. This is due to its sensitivity to young, massive stars and dust attenuation.

A.1. Fundamental plane: Massive spheroids with M_* > 1 × 10¹¹M_⊙

We find the simulated L-FPs are well defined, exhibiting relatively low scatter at all redshifts.

In Fig. A.1 we present the redshift evolution of L-FP slope coefficients b (upper panel) and c (middle panel), as well as the zero-point offset (Δa = a_z − a_zi; bottom panel), with comparisons to other simulations and observational studies. We find that for our massive galaxy sub-sample the evolution of the coefficient c represented by the solid red line agrees well with the Illustris-TNG results (dashed pink line; Lu et al. 2020). In contrast, the Horizon-AGN simulation (dashed green curve; Rosito et al. 2021) predicts larger absolute values of c, deviating from the other two simulations. Despite these differences in magnitudes, all simulations suggest near redshift-independence of the coefficient c. On the other hand, all three simulations predict roughly the same value of the coefficient b up to z = 2, with only minor difference within statistical uncertainty. An interesting finding is that they all predict a weak dependence of b on redshift in the sense that b monotonically decreases toward higher redshift, indicating that coefficient b is more sensitive to evolutionary effects than coefficient c. Further in bottom panel of Fig. A.1, we find that Δa (solid red curve) remains relatively flat consistent with zero. This trend suggests that the structural and dynamical scaling relations of massive galaxies were already largely in place by z ∼ 2.5, with little subsequent evolution in the FP zero-point.

It is important to note that ETG selection criteria vary across simulations and different observational studies, which may explain some of the quantitative discrepancies. For example, in Illustris-TNG simulation a galaxy is classified as an early type if its radial surface-brightness profile is better fit by a de Vaucouleurs law than by an exponential profile, and if its sSFR lies at least 1 dex below the star-forming main-sequence ridge.

Fig. A.1. Redshift evolution of the L-FP coefficients, b (top), c (middle), and the zero-point offset, Δa = az − azi (bottom), where the L-FP is defined as log Rhsm = a + blog σe + clog Ie. Results from this work are shown for two stellar mass-selected samples: a massive sub-sample with M* > 1 × 1011 M⊙ (solid red line) and the full sample with M* > 2 × 109 M⊙ (blue dash-dotted line). For comparison, we also show L-FP results from other cosmological simulations (HorizonAGN: dashed green line; TNG100: dashed magenta line) and observational studies (Jørgensen et al. 1996; van Dokkum & Franx 1996; Saglia et al. 2010; Hou & Wang 2015; Stockmann et al. 2021).

A.2. Fundamental plane: Slowly rotating spheroids with M_* > 2 × 10⁹M_⊙

We now turn to our analysis of a broader spheroidal galaxy sample with stellar masses M_* > 2 × 10⁹ M_⊙. Similar to the massive sub-sample, we find that the FP remains well-defined across all redshifts, from z = 0.625 up to z = 6, with relatively low scatter throughout. The evolution of the FP coefficients b and c for this lower-mass sample is shown in Fig. A.1 (blue curves). Compared to the massive galaxy sample, the lower-mass galaxies exhibit larger values of b and smaller absolute values of c, indicating mass dependence in the FP slopes. The coefficient c shows little evolution up to z ∼ 3 with an indication of rise up to z ∼ 6. In contrast, we do see a clear redshift evolution of the coefficient b for the broader galaxy sample. It decreases slowly as redshift increases beyond z ∼ 1.2 but now all the way to z ∼ 6. The value of b is about 40% larger for lower mass galaxies.

In the bottom panel of Fig. A.1, we show the redshift evolution of the zero-point offset (Δa) for galaxies with stellar masses M_* > 2 × 10⁹ M_⊙, alongside predictions from the Horizon-AGN (dashed green line) and TNG100 (dashed pink line) simulations for comparison. A clear increasing trend is observed in Δa for the lower-mass galaxy sample (blue line), particularly up to z ∼ 2.5 agreeing with Horizon-AGN results. However, there is weak decline in Δa beyond z > 2.5. This behaviour may reflect a transition in the dominant physical processes affecting galaxy structure at high redshift such as increased merger rates, evolving mass-to-light ratio, etc. Further investigation is needed to understand the physical mechanisms driving this trend.

All Tables

All Figures

	Fig. 1. Comparison of VCF for the various cosmological hydrodynamical simulations highlighting the importance of larger simulation box size. It can be seen that the VCF is significantly underestimated at zero separation (velocity dispersion), galaxy scale (r ≈ 1 h−1 Mpc) and cluster scale (r ≈ 10 h−1 Mpc) when the simulation size is less than 1cGpc
In the text

Fig. 2. Distribution of 3D half stellar mass radius (log Rhsm), SDSS (g − r) colour, stellar age in r band (log ager), and Sérsic index (nSérsic) versus stellar mass (log M*) for selected spheroidal galaxies (blue points), their massive sub-sample (${\mathrm{M}}_{\ast }>{10}^{11}{\mathrm{M}}_{\odot }$; orange points), and the full sample of galaxies (with ${\mathrm{M}}_{\ast }>2\times {10}^9{\mathrm{M}}_{\odot }$; black contours) from the HR5 simulations at z = 0.625, respectively. The dashed black line in the bottom right-hand panel indicates nSérsic = 1.5.

In the text

Fig. 3. Redshift evolution of SM-FP for the spheroidal galaxies in HR5. Note that galaxies are located along the diagonal line as the best-fit FP coefficients are used at each redshift. The relation moves to the right as galaxies grow in size with time. Left: For the massive sub-sample with M* > 1 × 1011 M⊙Right: For the sample with M* > 2 × 109 M⊙. To mitigate the overlap of points, black contours are plotted for the initial five redshifts, representing the central 90% of galaxies contained within.

In the text

Fig. 4. Redshift evolution of the SM-FP parameters: slope coefficients, b (top) and c (middle), and zero-point variation, Δa = az − azi (bottom). The results of the HR5 samples are indicated by a red curve for the massive sub-sample and a blue curve for the full sample. The errors on a, b, c, Δ, and Δ/Δr were estimated using bootstrap resampling for the massive sub-sample, while for the broader sample that includes lower-mass galaxies, they were determined with the sub-sampling method. Observational constraints from Hou & Wang (2016) and Bernardi et al. (2020) are shown as magenta and yellow points, respectively.

In the text

	Fig. 5. Redshift evolution of SM-FP for the spheroidal galaxies in HR5. The symbols are coloured according to the mass-weighted average stellar age (upper panels), Vrot/σe (middle panels), and sSFR = SFR/M*[yr−1] (bottom panel). The 1:1 line is depicted in black in each panel.
In the text

Fig. 6. Comparison between the L-FP and the SM-FP as a function of redshift. The panels show the evolution of FP coefficients b (top) and c (middle), and the zero-point offset, Δa = az − azi (bottom). Results are shown for two stellar mass-selected samples: a massive sub-sample with M* > 1 × 1011 M⊙ (red) and the full sample with M* > 2 × 109 M⊙ (blue). Solid lines correspond to the SM-FP, while dashed lines indicate the L-FP.

In the text

Fig. A.1. Redshift evolution of the L-FP coefficients, b (top), c (middle), and the zero-point offset, Δa = az − azi (bottom), where the L-FP is defined as log Rhsm = a + blog σe + clog Ie. Results from this work are shown for two stellar mass-selected samples: a massive sub-sample with M* > 1 × 1011 M⊙ (solid red line) and the full sample with M* > 2 × 109 M⊙ (blue dash-dotted line). For comparison, we also show L-FP results from other cosmological simulations (HorizonAGN: dashed green line; TNG100: dashed magenta line) and observational studies (Jørgensen et al. 1996; van Dokkum & Franx 1996; Saglia et al. 2010; Hou & Wang 2015; Stockmann et al. 2021).

In the text