© The Authors 2023

1. Introduction

Galaxies in cluster environments have different star formation histories and display different morphology types with respect to galaxies in the field. Generally speaking, galaxies in clusters are redder, more likely to be quenched, and exhibit an elliptical morphology (e.g., Dressler 1980; Butcher & Oemler 1984; Dressler & Oemler 1997; Peng et al. 2010; Lewis et al. 2002; Blanton et al. 2005; Alpaslan et al. 2015). Galaxy evolution is mainly affected by three physical processes related to the cluster environment. First, in a process known as “harassment”, tidal forces induced by the central halo (Gunn & Gott 1972; Toomre & Toomre 1972; Barnes & Hernquist 1992; Bournaud et al. 2004) or neighbouring flyby galaxies can strip away part of the stellar and gas particles of the satellite galaxies (Gunn & Gott 1972). Second, the interaction with the intracluster medium (ICM) known as “ram pressure” can strip away the hot gas and sometimes even the cold gas in the disk of the satellite galaxies (Gunn & Gott 1972; Abadi et al. 1999; Yun et al. 2019). Third, the accretion rate is much lower in the cluster environment and it could cause the cut-off of gas supply to a galaxy leading to “strangulation” (Larson et al. 1980; Balogh et al. 2000; Kawata & Mulchaey 2008). All three processes could ultimately lead to the cessation of star formation and affect the morphology evolution of the satellite galaxies in the cluster.

The timescale of the star-formation quenching in cluster environments and the relative contribution of the aforementioned physical processes are still under debate. The current gas content in cluster galaxies provides some direct clues. Maier et al. (2019b) calculated the fraction of star-forming galaxies of seven nearby clusters through Hα observation with Local Cluster Substructure Survey (LoCuSS). By comparing with the Millennium simulations, they suggested that star formation in cluster galaxies can last for 1–2 Gyr after their infall into the cluster, but then they quickly get quenched. This “slow-then-rapid quenching” is consistent with the strangulation scenario. On the other hand, Reynolds et al. (2022) found evidence that HI gas in the outer regions of the galaxies in the Hydra I cluster is partially removed. However, these galaxies still lie on the star-forming main sequence and gas removal is not yet affecting the inner star forming disks. In the same cluster, Wang et al. (2021) modeled the ram-pressure stripping strength through the HI mass from WALLABY observation. They found that the ram pressure stripping can significant change the total HI gas mass of satellite galaxies within 600 Myr after falling into the cluster.

Cosmological hydrodynamical simulations allow us to investigate the details of the infalling of galaxies into the cluster as well as of the gas removal and structure formation of galaxies. Joshi et al. (2020) found that the pre-exisiting stellar disks in cluster galaxies are largely disrupted by impulsive tidal shocking and stripping at pericentres within 0.5–4 Gyr after falling into the massive clusters with M₂₀₀ ∼ 10^{14 − 14.3} M_⊙ in the Illustris TNG100 simulation. Meanwhile, all the cluster galaxies quenched their star formation after the disk disruption. The quenching timescale could be 1–1.5 Gyr for gas-poor massive galaxies, while the star formation could last for ∼4.5 Gyr for gas-rich galaxies and a dynamically cold disk could regrow. Similar results apply also to Fornax-mass analogues in the TNG50 simulations (Galán-de Anta et al. 2022)

The detailed study of S0 galaxies in TNG100 simulation shows two major formation scenarios (Deeley et al. 2021): S0 galaxies in clusters formed via starvation or stripping, whereas S0s in field formed via mergers. This finding is also supported by observations (Coccato et al. 2020, 2022). Star formation quenching is different in these two formation scenarios (Deeley et al. 2021): when a blue and gas-rich galaxy falls into a cluster, its gas in the outer regions is quickly stripped, while the star formation in the inner regions could last for a long time, whereas in case of galaxy–galaxy merger, galaxies start quenching from the inner regions, and star formation continues in the outer ring-like regions.

These studies suggest that cluster environment may significantly affect the internal dynamical structure of galaxies, especially the persistent growth of dynamically cold disks. The fraction of stars in cold disks and their age distributions may imprint the long-term star-formation quenching process happening in cluster galaxies. We expect that the disk stars are older in the inner regions and younger in the outer ones for a normal inside-out growth (e.g., Bird et al. 2013; Stinson et al. 2013). Cluster quenching may turn over the age gradient in the stellar disks, the age difference between inner and outer regions can thus tell us the star formation duration in galaxies after falling into the cluster. Although some other studies show that complicated physical processes in field galaxies might cause similar age gradient (Lagos et al. 2022).

In the past decades, integral-field unit (IFU) surveys, such as SAURON (de Zeeuw et al. 2002), ATLAS^3D (Cappellari et al. 2011), CALIFA (Sánchez et al. 2012), SAMI (Bryant et al. 2015), and MaNGA (Bundy et al. 2015), have provided us the stellar kinematics, age, and metallicity maps for thousands of nearby galaxies. Early-type galaxies (ETGs) show a large variety of kinematical structures (Cappellari et al. 2007; Emsellem et al. 2007, 2011), some lenticular galaxies have a rapidly rotating disk similar to that of spiral galaxies. This suggests that gas removal by cluster environment may play a major role in transforming spirals into lenticulars (Cappellari 2016). Overall, ETGs are found to be old and with shallow age gradients (McDermid et al. 2015; González Delgado et al. 2015; Zibetti et al. 2020), while late-type galaxies (LTGs) generally have negative age gradients, which is consistent with the inside-out scenario (González Delgado et al. 2015). The environmental dependence of the galaxy age gradient is still controversial. A positive age gradient was found in SAURON and MaNGA samples (e.g., Goddard et al. 2017; Li et al. 2018; Kuntschner et al. 2010), with a mixture of galaxies in cluster and field environments. On the other hand, there is no significant difference found for age gradient of ETGs in different environments (Zheng et al. 2017; Ferreras et al. 2019; Santucci et al. 2020), with the above surveys covering the inner 1R_e ∼ 2R_e for most galaxies.

The Fornax3D survey (Sarzi et al. 2018) measured high-quality age and metallicity maps covering at least the inner 2R_e of 23 ETGs in the Fornax cluster. These galaxies show significant positive age gradients between R_e and 2R_e (Spavone et al. 2022), while they have negative metallicity gradients, which are shallower than the control sample galaxies in the field. On the other hand, positive age gradients together with negative metallicity gradients are widely found in the dwarf galaxies of the Local Group (e.g., Koleva et al. 2011; Kirby et al. 2011; Zhuang et al. 2021). These positive age gradients could be the direct consequence of ram pressure stripping (Genina et al. 2019; Deeley et al. 2021), although others argue theses are the result of gas-rich mergers for dwarf galaxies (Yozin & Bekki 2012).

The observation of external galaxies provides integrated properties along the line-of-sight. The stellar kinematics and populations of a galaxy are a combination of different dynamical components, displaying a range of possible physical origins. As learned from the simulations aforementioned, the luminosity fraction and age distribution of stellar disk might be a good probe of environmental affects. However, there are relatively small disk fractions in ETGs (Zhu et al. 2018a). The age and metallicity maps of the whole galaxy may not reveal the properties of the disk component.

Dynamical models are powerful tools to uncover the galaxies’ underlying mass profiles, as well as the internal 3D structures which could lead to physical-motivated decomposition of different components. There are a few methods widely used for modelling the IFU kinematic data of external galaxies, including the Jeans-Anisotripic-MGE (JAM) models (e.g., Cappellari 2008; Li et al. 2017), the particle-based Made-to-Measure (M2M) models (e.g., Syer & Tremaine 1996; de Lorenzi et al. 2007; Long & Mao 2010, 2012; Zhu et al. 2014), and the orbit-superposition Schwarzschild method used in this work (Schwarzschild 1979). The Schwarzschild’s orbit superposition method can be used in different geometries such as the spherical systems (Richstone & Tremaine 1985; Breddels et al. 2013; Kowalczyk et al. 2017), axisymmetric systems (Cretton & van den Bosch 1999; Gebhardt et al. 2000; Valluri et al. 2004; Cappellari et al. 2006a; Thomas et al. 2007; Saglia et al. 2016), and triaxial systems (van den Bosch et al. 2008; Neureiter et al. 2021). The orbit-superposition model reconstructs the backbone of galaxies without ad-hoc assumptions of the underlying distribution functions, it has been widely used to uncover galaxies’ underlying dark matter distributions (e.g., Cappellari et al. 2006b; Yang et al. 2020), central black hole mass (e.g., van der Marel et al. 1998; Cretton & van den Bosch 1999; Verolme et al. 2002; Gebhardt et al. 2003; Valluri et al. 2004; Krajnović et al. 2009; Ahn et al. 2018; Thater et al. 2022b), and internal stellar orbit distributions (e.g., Zhu et al. 2018a,b; Jin et al. 2020). Recently it has been modified to include the barred structures explicitly, thus also for uncovering the bar pattern speed (Vasiliev & Valluri 2020; Tahmasebzadeh et al. 2022).

Based on the Schwarzschild method, a population-orbit superposition method (Poci et al. 2019; Zhu et al. 2020) was recently developed by tagging age and metallicity to the orbits, we can thus simultaneously fit the stellar kinematic, age, and metallicity maps from IFU survey, and obtain the internal stellar orbit distribution of galaxies as well as the age and metallicity distributions. This allows for a physically motivated chemo-dynamical decomposition of galaxies. In this work, we apply the population-orbit superposition method to galaxies in Fornax cluster with the IFU data obtained with MUSE/VLT in the context of the Fornax3D survey (Iodice et al. 2019a). Using this modeling technique, age-dispersion profiles of dynamically decomposed cold disks in a few edge-on galaxies have been obtained (Poci et al. 2019, 2021) and a dynamically decomposed “hot inner stellar halo” has been used to weigh and time the ancient massive mergers the galaxies experienced in NGC 1380 and NGC 1427 (Zhu et al. 2022).

In this paper, we set our focus on dynamically cold disks. By studying the luminosity fraction, age, and metallicity of these stars, we aim to quantify how the cluster environment has affected the formation of cold disks. This results will also lead to a direct and in-depth comparison with galaxies in cosmological simulations and, thus, to a better understanding of the physical processes driving galaxy evolution in cluster environments. This work is undertaken in the context of ΛCDM cosmology, with density parameters of Ω_m = 0.3089, Ω_Λ = 0.6911, Ω_b = 0.0486, normalization σ₈ = 0.8159 and spectral index n_s = 0.9667 (Planck Collaboration XVI 2014).

The paper is organized as follows. We introduce the data set in Sect. 2. We describe the population-orbit superposition model and its relevance to the orbital decomposition in Sect. 3. We show the orbital decomposition for all galaxies and present the surface-brightness, age, and metallicity radial profiles of each component in Sect. 4. We present an explanation of how the cold disk age can be used as a novel proxy for galaxy infall time into a cluster in Sect. 5. We further show the dependence of cold disk properties on galaxy stellar mass and cluster environment in Sect. 6. We discuss the results in Sect. 7 and present our summary in Sect. 8.

2. Sample and data

We set our focus on galaxies in the Fornax cluster, which has a virial radius of R_vir ∼ 0.7 Mpc, virial mass of M_vir ∼ 7 × 10¹³ M_⊙, and a distance of D ∼ 20 Mpc (Diaferio 1999; Drinkwater et al. 2001).

We used the deep photometric data of the Fornax Deep Survey (FDS; Venhola et al. 2018). FDS observed all the galaxies in the area of 9 square degrees around the core of the cluster, with the VLT Survey Telescope down to a surface-brightness level of 27 mag arcsec⁻² in the r band (Iodice et al. 2016).

We used the spectroscopic data of the Fornax3D survey (Sarzi et al. 2018). Using the MUSE instrument on the VLT, Fornax3D observed all the 23 ETGs and 10 LTGs within the virial radius of the cluster down to 25 mag arcsec⁻² in B band. MUSE has a field-of-view of 1 × 1 arcmin² and spatial scale of 0.2 arcsec pixel⁻¹. The wavelength covers the range of 4650 Å–9300 Å, with a resolution of FWHM_7000&Å = 2.5 Å. Twenty-eight galaxies were observed by more than one MUSE pointing. The data extend to at least 2R_e for 20 galaxies.

2.1. Stellar kinematics maps

To obtain reliable stellar kinematics, we first bin the spectra in nearby pixels to reach a target signal-to-noise ratio (S/N) taken a Voronoi binning scheme (Cappellari & Copin 2003). The stellar kinematics of all the galaxies in the Fornax3D survey was published in Iodice et al. (2019a) with target S/N = 40. The only difference in the kinematic maps we use here is that different S/N for different galaxies were chosen as indicated in Table 1.

We took a target S/N = 100 or 200 for the bright galaxies with more than one MUSE pointing and S/N = 60 or 40 for the remaining galaxies. This results in a number of bins from ∼100 to ∼1000 for each galaxy with a spatial resolution from ∼100 pc in the inner region to ∼1 kpc in the outer faint regions. This choice is convenient for dynamical modeling without too many bins, while the spatial resolution is still good enough to meet our science goals for resolving sub-kpc scale structures.

The stellar kinematics was then extracted by applying the pPXF full-spectral fitting (Cappellari & Emsellem 2004) to the wavelength range 4750–5500 Å. This yielded high-quality maps of the stellar mean velocity, V, velocity dispersion, σ, and higher order velocity moments parameterized through the Gauss-Hermite (GH) coefficients h₃ and h₄ (Gerhard 1993; Rix et al. 1997).

Such stellar kinematic maps were derived in a consistent way for the 15 ETGs and 5 LTGs included in this paper. We did not consider the 8 ETGs – FCC 213 (i.e., the BCG), FCC 090, FCC 184, FCC 190, FCC 193, FCC 219, FCC 277, and FCC 310 – due to existence of bright foreground stars, limited data coverage or the presence of strong bar, and the 5 LTGs FCC 113, FCC 176, FCC 267, FCC 285, and FCC 306 due to the presence of a strong bar or irregular morphology. We excluded the strongly barred galaxies because our current dynamical models cannot fit their kinematic maps well with bars that are not included explicitly. However, we did keep the three weakly barred galaxies (FCC 179, FCC 182, and FCC 263), since their stellar kinematic maps are not strongly affected by the bar (see figures in Appendix C) and our current kinematic models can still provide reasonably good fits. Some basic properties of the final sample of galaxies are listed in Table 1.

2.2. Age and metallicity maps

We used the age and metallicity maps derived in Martín-Navarro et al. (2021). Here, we briefly describe how they were derived.

To obtain high-quality maps of the age and metallicity, the spectra were spatially rebinned to reach a minimum of S/N = 100 for each galaxy. We note that the spaxels with S/N < 5 were not included in this process in order to retain the highest quality of data. Thus, for most of the galaxies, only the inner region map was available. The single stellar population synthesis models (SSP) by Vazdekis et al. (2016) based on the MILES stellar library by Sánchez-Blázquez et al. (2006) at a constant spectral resolution of 2.51 Å (FWHM; Falcón-Barroso et al. 2011) were used for the analysis.

The luminosity-weighted age and metallicity maps were derived in a two-step spectral fitting of the wavelength range between 4800 Å and 6400 Å, considering the initial mass function (IMF) variation. First pPXF was run with an un-regularized combination of MILES SSPs in order to measure V and σ. A second pPXF was run again with stellar kinematics fixed to that measured in the first step, and this time regularizing the age–metallicity–IMF slope plane. The mean luminosity-weighted age and metallicity of the spectra were measured using this regularized second pPXF run. The metallicity, Z, is added up linearly rather than in the commonly used logarithmic way. Templates with a variable IMF were intentionally included to account for the possible effect of the IMF in the age determination. The abundance ratio [α/Fe] was not fitted, but we used the so-called base MILES models which inherit the [α/Fe]–[M/H] relation of the solar neighborhood (Vazdekis et al. 2015). This choice avoids non-local equilibrium uncertainties introduced by the theoretical response functions needed to compute stellar population models with variable abundance ratios, uncertainties that can be particularly problematic for Balmer lines and thus for the age determination. We note that in Martín-Navarro et al. (2021), the metallicity maps were also derived in a third step with index fitting when determining the IMF. However, we used the luminosity-weighted metallicity maps derived from the pPXF fitting, in a consistent way as the luminosity-weighted age maps.

Alternatively, we have another version of the age and metallicity maps for FCC 153, FCC 167, FCC 170, and FCC 177, which was obtained by pPXF fitting regularizing age–metallicity–[α/Fe] assuming a Kroupa IMF (Pinna et al. 2019b,a).

We compared the two versions of age and metallicity maps for these four galaxies, they generally show the same age and metallicity gradients but with a systematic offset of 1–2 Gyr in age. This offset will not affect our results. The version of the age and metallicity maps by Martín-Navarro et al. (2021) is used consistently for all galaxies in this paper.

Finally, we have the age and metallicity maps for 16 galaxies, including all the 15 ETGs and the LTG FCC 179 only. The remaining 4 LTGs of our sample appear to be dominated by young stars at all radii, for which the uncertainties of age and metallicity derived from the spectra fitting are large (Ge et al. 2018), we will build their dynamical models without associating age and metallicity to the stellar orbits.

3. Population-orbital superposition method

We follow Zhu et al. (2020)¹ in constructing population-orbit superposition models that simultaneously fit the luminosity density, kinematic, age, and metallicity maps of each galaxy. From the best-fit models, we obtain the internal stellar orbit distribution associated with the age and metallicity distributions. Based on the model, we then decomposed each galaxy into a dynamically cold disk and a dynamically hot non-disk component, based on circularity distribution of the orbits, and extracted the face-on surface-brightness, age, and metallicity radial profiles for each component separately. We used FCC 177 as an example to illustrate the whole process as shown in Figs. 1 and 2.

Fig. 1. Best-fit population-orbit superposition model of FCC 177. From top to bottom: maps of the data, model, and residuals (model−data). From left to right: maps of the surface-brightness, SB, mean velocity, V, velocity dispersion, σ, Gauss-Hermite coefficients h3 and h4, light-weighted age and metallicity of the stars. Similar plots for the other sample galaxies are shown in Figs. B.1–B.19.

Fig. 2. Orbital decomposition of FCC 177. Probability density distribution p(r, λz) in the left panel, age distribution p(r, t) in the central panel, and metallicity distribution p(r, Z) in the right panel of the stellar orbits in the phase space of time-averaged radius, r, vs. circularity, λz. The probability densities are normalized to unity within the data coverage. All the distributions are averages of multiple best-fit models that fall within the 1σ confidence level of the model hyper-parameters. The dashed line marks our orbit-based division into two components: a dynamically cold disk component (λz ≥ 0.8) and a dynamically hot non-disk component (λz < 0.8). The shadow regions are beyond the data coverage of age and metallicity maps. Similar plots for the other galaxies are shown in Figs. B.1–B.19.

The construction of a population-orbit superposition model consists of four steps. First, constructing the model of gravitational potential. Second, calculating the orbital library under the gravitational potential. Third, fitting the luminosity density and kinematic maps by weighting the orbits. Fourth, fitting the age and metallicity maps by “coloring” the orbits. The method has been presented in great detail and carefully validated in Zhu et al. (2020). Here, we just briefly describe some of the key steps.

3.1. Constructing the gravitational potential

The gravitational potential combines the contribution of the stars, dark matter halo, and a central black hole.

To construct the stellar mass distribution, we use the multi-Gaussian expansion (MGE; Cappellari 2002) to fit the r-band image from FDS. The galaxies in our sample have a variety of morphological and kinematic properties, which can be classified in three categories: galaxies which are generally flat and with strong rotation indicating an extended disk, galaxies which are generally round but flatter, and with stronger rotation in the inner regions indicating an inner disk, galaxies with no clear evidence of a disk. For the first two categories, we consider an axisymmetric solution by adopting a constant positional angle (PA) for the Gaussians to match that of the photometric and kinematic signature of the inner or outer disk component. Photometric PA from the disk region is taken, which is generally consistent with the kinematic PA in that region. For the third category, we still consider axisymmetric solution, but we did not put any extra constrain on the PA; thus, the photometric PA directly derived from the whole galaxy image was taken. The best-fit Gaussian parameters are listed in the appendix for all the galaxies.

Then we de-project the 2D Gaussian components to 3D by assuming a set of viewing angles (θ, ψ, ϕ), where θ and ϕ define the orientation of the line of sight with respect to the principal axes of the galaxy and ψ is chosen to specify the rotation of the galaxy around the line-of-sight in the projected sky-plane. By combining the 3D Gaussian components, we obtained the 3D luminosity density distribution of the galaxy and we multiply a stellar mass-to-light ratio M_*/L to get the 3D mass density distribution. The gravitational potential could then be calculated by the classical Chandrasekhar formula (van den Bosch et al. 2008).

In practice, we do not directly use the three viewing angles as free parameters. In van den Bosch et al. (2008), three intrinsic shape parameters (p, q, u) are used instead of the viewing angles (θ, ψ, ϕ), where q = Z/X, p = Y/X, and u = X′/X, where X, Y, Z are the intrinsic long, intermediate, and short axis of the galaxy and X′ is the projected major axis. The conversion of the two sets of parameters follows Eq. (10) in van den Bosch et al. (2008). Exploring the space of intrinsic shape parameters is more efficient than that of the three viewing angles. For instance, the deprojection of an axisymmetric system would have no constrain on the parameter ϕ but have a finite axis-ratio between Y and X. In this work, we follow the approach by van den Bosch et al. (2008) and allow for some degree of triaxiality of the galaxy by setting a non-unity u.

We use a spherical Navarro–Frenk–White (NFW; Navarro et al. 1996) dark matter halo, with one free parameter dark matter virial mass M₂₀₀, while concentration C is fixed by the correlation between M₂₀₀ and C (Dutton & Macciò 2014). Although there is a large scatter in the M₂₀₀ − C plane for real galaxies, the choice of a fixed C will not significantly affect our results because the two parameters are degenerated and will not be constrained separately with our kinematic data covering out to 2 − 3R_e. The potential also includes a central black hole characterized by a Plummer potential (van den Bosch et al. 2008), with the back hole mass, M_•, fixed according the M_• − σ relation (Kormendy & Ho 2013). This choice does not affect our results, because the black hole sphere of influence is mostly not resolved by our kinematic data and we cannot directly constrain the black hole mass.

In summary, we have five free “hyper-parameters,” for the gravitational potential: the mass-to-light ratio, M_*/L, three parameters on the intrinsic shape of stellar distribution, p, q, and u, and the dark matter virial mass, M₂₀₀.

3.2. Computing the orbit library

For each model, with a set of hyper-parameters, we calculated an orbit library with tens of thousands orbits. The orbits sampling follows the way described in van den Bosch et al. (2008).

We first sampled regular orbits according to a separable triaxial potential. The orbits are sampled from the three integrals of motion: energy, E, second integral of motion, I₂, and third integral of motion, I₃ (Binney & Tremaine 2008). We sample two sets of 55 × 11 × 11 combinations of (E, I₂, I₃) as initial conditions, which include co-rotating and counter-rotating orbits.

Box orbits are crucial for supporting the triaxial structures, we sampled another set of box orbits by constructing the initial conditions on equipotential surfaces with the energy, E, two spherical angles, θ and ϕ, which gives another set of 55 × 11 × 11 orbits.

The number of orbits, especially across E, we sample here is larger than that used in previous works (e.g., van den Bosch et al. 2008; Zhu et al. 2018c; Jin et al. 2019). For instance, 21 × 7 × 7 orbits were used in fitting the CALIFA data (Zhu et al. 2018c). Because the kinematic data used in this work have larger spatial coverage and higher spatial resolution than those from the previous IFU surveys, more freedom is required by the model to fit the data. To reduce the Poisson noise of the model, we dither every orbit by slightly perturbing the initial conditions to give 5³ orbits.

3.3. Fitting stellar luminosity density and kinematics

We fit the luminosity density and stellar kinematics of the galaxy by weighting the orbits. The 2D surface-brightness, 3D luminosity density deprojected from the MGE, and kinematic maps are used as model constrains. The kinematic maps include the line-of-sight mean velocity, V, velocity dispersion, σ, Gauss-Hermite coefficients, h₃ and h₄. It is worth noting that we did not fit V and σ maps directly, but the Gauss-Hermite coefficients h₁, h₂, h₃, and h₄. We extracted similar luminosity and kinematic maps from the model superposed by orbits, then obtained a solution of the orbit weights by minimizing the χ² between the data and model using a non-negative least squares (NNLS) method (Lawson & Hanson 1974, 1995).

We used an optimized grid searching process (Zhu et al. 2018b) to adjust the free hyper-parameters of the gravitational potential. We started with a model with an initial guess of the hyper-parameters, then we performed an iterative searching process with intervals of 0.05, 0.02, 0.01, 0.01, and 0.05 for the hyper-parameters, M_*/L, p, q, u, and log M₂₀₀. We note that the concentration parameter, C, in the NFW dark matter halo is fixed by the M₂₀₀ − C relation. After the previous sampled models were completed, we selected the best-fit models with and sampled the new models around the selected ones. We continued the iterative process until an area of χ² minimum was found and all models within 3-σ confidence level around the minimum χ² were calculated. At the end, we calculated a few hundred models for each galaxy, the resulting parameter grid for FCC 177 is shown in Fig. C.1. The models with least-χ² are selected as the best-fit model. In Fig. 1, we show the best-fit model of FCC 177, where the model matches the observed kinematic data in great detail. Considering the model uncertainties, we defined the 1σ confidence level in a similar way as Zhu et al. (2018a):

(1)

where n_GH = 4 is the number of stellar kinematic moments and N_obs is the number of bins in the kinematic data. We adopt f = 1 for the galaxies with N_obs < 200, similar to the previous models for the CALIFA data; for galaxies with a much larger value for N_obs, we found that the χ² fluctuation caused by numerical noise is much higher and thus we adopted f = 4, following the results of a bootstrapping analysis².

Once we fit the stellar kinematics, we obtained the intrinsic stellar orbit distribution of the galaxy. We use the circularity, λ_z, and time-averaged radius, r, to characterize the orbits, where λ_z is defined as the orbital angular momentum around the z direction, normalized by the maximum that is allowed by a circular orbit with the same binding energy. Here, λ_z ∼ 1 represents highly rotating short-axis tube orbits and the λ_z ∼ 0 represents mostly long-axis tube or box orbits. The stellar orbit distribution from the best-fit model of FCC 177 is shown in the left panel of Fig. 2.

3.4. Tagging orbits with age and metallicity

Next, we fit the age and metallicity maps by tagging age and metallicities to the orbits. We took the models within the 1σ confidence level selected by the kinematic fitting as given by Eq. (1). For each model, we have the stellar orbit distribution. We applied a Voronoi binning scheme to the orbits in the phase-space of λ_z vs. r and decompose them into ∼100 orbital bundles. Orbits with similar r and λ_z are included in the same orbital bundle and each bundle has a minimum orbital weight of 0.005.

After the Voronoi binning, we assume that each orbital bundle k has a single value of age, t_k, and metallicity, Z_k. Then the age, , and metallicity, , of each observational aperture, i, can be expressed as:

(2)

(3)

where N_bundle is the total number of orbital bundles and is the luminosity contribution of orbital bundle, k, at an aperture, i.

We applied a Bayesian statistical analysis (Python package pymc3; Salvatier et al. 2016) to first solve t_k by fitting the observed age map and then to solve Z_k by fitting the observed metallicity map. To use the Bayesian theorem to compute the posterior probability of a model, it is necessary to include the prior probability and data likelihood. For the prior probability of t_k, we adopted a bounded normal distribution, where the mean value μ(t_k) is randomly sampled around the average of the observed . We adopted a student’s t-distribution for the data likelihood (Salvatier et al. 2016), which allows for some outliers in the data and results in a robust fitting. For each t_k, we ran a chain with 2000 steps and take the last 500 steps to calculate the mean and 1σ values of t_k.

For a few galaxies, we did not get a good match of the age maps with the above fitting process. Therefore, we tried a second fitting round with μ(t_k) in the bounded normal distribution chosen as the t_k value obtained from the first fitting round (Zhu et al. 2020). We obtained a good match of age maps for all the galaxies after the second round of fitting.

Then we fit the metallicity map. For the prior probability of Z_k, we adopt a log-normal distribution, where the mean value μ(log(Z_k)) is set by a build-in age-metallicity relation using t_k obtained from the fitting of age map and following Zhu et al. (2020). The student’s t-distribution is adopted again for the data likelihood, and we run a chain with 2000 steps and take the last 500 steps to calculate the mean and 1σ values of Z_k.

The best-fit model of FCC 177 is shown in Fig. 1. At the end, we obtained models that simultaneously fit well the surface-brightness, stellar kinematics, age, and metallicity maps of all the galaxies. We note that the age and metallicity maps we used have different binning schemes and have less data coverage than the kinematic data. For this reason, the age and metallicity of the orbits beyond the data coverage are not well constrained.

3.5. Orbital decomposition

In Fig. 2, we show the probability density, age, and metallicity distributions of the orbits of the best-fit model of FCC 177 in the phase-space of λ_z vs. r. We consider only the age and metallicity distributions within the data coverage (r ≲ 45 arcsec) to be reliable. The models are usually more noisy (than smooth) in the phase-space of λ_z vs. r. As we tested with mock data in Zhu et al. (2020), not all the substructures in the model are real, but we can trust the model considering the distributions of orbits and of these age and metallicity in a statistical way.

Based on the stellar orbit distribution, we coarsely decompose the galaxy into two components: the orbits with λ_z ≥ 0.8 are taken as a dynamically-cold component which is usually a disk. All the remaining orbits with λ_z < 0.8 are considered as a dynamically-hot non-disk component. We note that here the dynamically cold component is defined as in Zhu et al. (2018c), whereas the dynamically hot component actually includes their dynamically warm, hot, and counter-rotating orbits. In Fig. 3, we show the radial profiles of surface-brightness, age and metallicity of the whole galaxy for FCC 177, as well as those profiles for the dynamically cold and hot components, separately. At r ∼ 0.3R_e, the cold disk is about 2 Gyr younger and significantly more metal-rich than the hot component and the age difference has risen to ∼4 Gyr at r ∼ 0.7R_e. However, there is a young and metal-rich nuclear star cluster in the galactic center (r ≲ 0.2R_e) which contributes to the hot component. The luminosity-weighted mean age and metallicity of the hot component are dominated by this young and metal-rich central regions, as a result, the difference of mean age (metallicity) of the cold disk and the hot component becomes smaller than that revealed in the profiles at r ≳ 0.3R_e. Within R_e, the stars on the cold disk orbits have an luminosity-weighted average age of 6.7 Gyr and an luminosity-weighted metallicity of 1.25 Z_⊙ while the stars on dynamically warmer orbits have an average age of 8.1 Gyr and an average metallicity of 0.81 Z_⊙.

Fig. 3. Surface-brightness, age, and metallicity radial profiles of the whole galaxy, the dynamically cold disk, and the dynamically hot non-disk component for FCC 177. We show the profiles of the whole galaxy (left panels) with the blue solid curve, and the profiles of dynamically cold disk and dynamically hot non-disk component (right panels) with the blue solid curve and black dashed curve. The shadowed areas indicate the scatter of the profiles of models that fall within the 1σ confidence level. The age and metallicity profiles of the dynamically cold disks are considered reliable and shown in the regions where the dynamically cold disk contributes at least 10% of the total surface brightness.

4. Results of the orbital decomposition for the sample galaxies

We have 20 galaxies with orbit-superposition models. For 16 of those (15 ETGs and one LTG, FCC 179), we also have age and metallicity maps and subsequently build a population-orbit model as well. In this section, we show the surface-brightness, age, and metallicity radial profiles of the dynamically cold disk component and dynamically hot non-disk component for all 16 galaxies with their orbits tagged with age and metallicity. For the four LTGs without orbits tagged with age and metallicity, we are limited to the surface-brightness profiles of each component.

4.1. Surface-brightness, age, and metallicity radial profiles

As illustrated by the best-fit model of FCC 177, we decomposed each model into a dynamically cold disk (with λ_z ≥ 0.8) and a dynamically hot non-disk component (including all the remaining orbits with λ_z < 0.8). Then, we reconstruct the 3D luminosity, age and metallicity distributions of each orbital component. We project each orbital component to be face-on and then extract the radial profiles of surface-brightness, age, and metallicity from the projected maps.

We note that the model already takes into the account the information for the 2D maps of the observational data. When we extract the radial profiles from the face-on projection, the radial profile along different directions are the same for the disk which is axisymmetric, while very similar for the hot non-disk component which could be moderately triaxial. For each galaxy, we have tens to hundreds of models within the 1σ confidence level. The aforementioned radial profiles are extracted for each model. We take the average of these profiles from models within the 1σ confidence level as the mean profile, and their scatter as the 1σ uncertainty.

There are significant variations of internal structure from galaxy to galaxy. We define the luminosity fraction of cold disk as:

(4)

where f_k is the luminosity fraction of orbital bundle k. The cold-disk fraction of a galaxy varies with radius, we calculated f_cold(r < R_e) and f_cold(r < 2R_e) with orbits within R_e and 2R_e, respectively.

We classify the galaxies into two categories: 11 of them have a relatively large cold-disk fraction and their cold disks have an extended surface-brightness, whereas the other 9 galaxies have a relatively small cold-disk fraction and their cold disks are concentrated in the very inner regions. Although the two categories are classified by eye, this is fairly consistent with a separation based on the cold-disk fraction f_cold(r < 2R_e): the galaxies in the first category generally have f_cold(r < 2R_e) > 0.1, most galaxies in the second category have f_cold(r < 2R_e) < 0.1, but a few > 0.1. We take this visual classification for practical reasons that we can obtain reliable age and metallicity gradients for the extended cold disks in the following analysis, but not for the concentrated ones.

The age and metallicity observed at any radius is a combination of different components. Our model tests in Zhu et al. (2020) show that we cannot recover the age and metallicity of a component at radii where it contributes ≤10% to the total luminosity. For the galaxies with extended disks, the surface-brightness of the cold disk component contributes at least 10% of the total surface-brightness over a radial region extending to at least one R_e. Reliable age and metallicity profiles of the disk are thus obtained over extended radial regions. On the contrary, for the galaxies with small cold-disk fractions and concentrated cold components, the age and metallicity radial profiles of the cold component are only derived in inner radial regions.

In Fig. 4, we show the radial profiles of surface-brightness of the whole galaxy and of the dynamically cold and hot components, for all the 20 galaxies, as well as the age and metallicity radial profiles for the 16 galaxies with their orbits tagged with age and metallicity. The galaxies with extended cold disks are shown in the left panels, while those with concentrated cold components are shown in the right panels. The dynamically cold and hot components are well separated with different surface-brightness, age, and metallicity radial profiles, which will be discussed in the following sections.

Fig. 4. Surface-brightness, age, and metallicity radial profiles of the galaxy, dynamically cold disk, and dynamically hot non-disk component for the sample galaxies. We divided the galaxies into two groups: 11 galaxies with a dynamically cold disk extended out to large radii (left panels) and 9 galaxies with a dynamically cold disk concentrated in the inner regions (right panels). The radial profiles are color-coded by the galaxy stellar mass. The shadowed areas indicates the scatter of the profiles of models that fall within the 1σ confidence level of the model hyper-parameters. Dashed lines refer to the four LTGs without orbits tagged with age and metallicity. The age and metallicity radial profiles of the dynamically cold disks are considered reliable and shown in the regions where the dynamically cold disk contributes at least 10% of galaxy surface-brightness. Therefore, these profiles are available only for a limited radial range for the galaxies with a concentrated cold disk.

We note that the surface-brightness profiles of the cold disks do not follow the traditional exponential law: the surface-brightness in the central regions is much lower than the inward extrapolation of the exponential profile. This finding is consistent with the results of orbital decomposition of the CALIFA galaxies (Zhu et al. 2018a) and with those from direct spectrum fitting, where bulge and disk are separated through their stellar populations (Breda et al. 2020).

4.2. Quantitative description of the radial profiles

We quantify the luminosity fraction f_cold of the cold disk component within R_e (or 2R_e), with respect to the cumulative luminosity of the galaxy within that radius. The 1σ uncertainty of f_cold is calculated by the scatter of the models within the 1σ confidence level. We calculate f_cold for all the 20 sample galaxies.

For the 16 galaxies with their orbits tagged with age and metallicity, we obtain the mean value and gradient of the radial profiles of age and metallicity of each component. We calculated the mean age and mean metallicity of the cold disk by

(5)

(6)

where t_k is the age of the orbital bundle, k. The lumonisity-weighted mean age ⟨t_hot⟩ and metallicity ⟨Z_hot⟩ of the hot component with the orbits λ_z < 0.8 are calculated in a similar manner. Note that the mean age is luminosity-weighted. Due to the different surface-brightness profiles of the cold and hot component, the difference in mean age of these two components sometimes appears differently from that seen in their age profiles, similar to the case of FCC 177 (as we explained in Sect. 3.5).

The gradients of age and metallicity are calculated as follows. We first uniformly interpolated 100 data points along the radial profile within the data coverage to smooth the profile, then we calculated the linear slope (gradient) between adjacent data points. In the end, we took the average of the gradients from the data points. We obtain the uncertainty of the gradient by bootstrapping within the shadowed regions of the profiles as shown in Fig. 4. All these parameters, together with other basic parameters obtained from the population-orbit superposition models, are included in Table 2.

5. Cold-disk age as a proxy of galaxy infall time into a cluster

5.1. Calibration from simulations

It is hard to accurately estimate the infall time into the cluster for each galaxy from observations, although we can statistically infer the likelihood of being ancient or recent infallers from their projected position and line-of-sight velocity in the cluster (Iodice et al. 2019a). By analyzing the cluster galaxies in the cosmological simulation Illustris TNG50 (Pillepich et al. 2019; Nelson et al. 2019), we find that the cold-disk age is tightly correlated with the infall time of the galaxy into the cluster, as a result of star formation quenching in disks associated with galaxy infall into the cluster (Ding et al., in prep.).

To make direct comparison between the Fornax and TNG50 galaxies. We defined cold disks of TNG50 galaxies in exactly the same way as here for the Fornax cluster galaxies, with the probability density distribution of stars in the phase-space of circularity λ_z versus r calculated from the 6D position-velocity information of particles in the simulation. We note that we use λ_z and r of the particles’ orbits, not directly the instantaneously values of each particle, following previous work in comparing the stellar orbit distribution of observed and simulated galaxies (Zhu et al. 2018c). The cold-disk fraction and cold-disk age are thus defined exactly the same way as described in Sect. 4. All galaxies in the 14 clusters with virial mass M₂₀₀ > 10^13.3 M_⊙ in the TNG50 simulations are chosen. We define the galaxies’ infall time into the cluster as the time when it first reaches the virial radius r₂₀₀ of the cluster at that time. For the pre-processed galaxies, we define their infall time as their first time of falling into the pre-cluster.

We found that the correlation between the cold-disk age and the infall time of the galaxy into the cluster does not strongly depend on the cluster mass. We divided all the galaxies in the 14 clusters into four mass bins with log₁₀M_*/M_⊙ ∈ (8, 8.6), (8.6, 9.2), (9.2, 10), and (10, 11), respectively. In the left panel of Fig. 5, we show the running median and ±1σ scatter of the correlation in the four mass bins color-coded by the stellar mass. The 1σ scatter is 0.63, 0.63, 0.94, and 1.34 Gyr from the low to high mass bins. The correlation is tighter for the least-massive galaxies. The scatter becomes significantly larger in the most massive galaxies because a significant fraction of them are already quenched or had little gas left before they fell into the cluster.

Fig. 5. Infall time into the Fornax cluster of the sample galaxies for which we obtained the age of the dynamically cold disk. Left panel: correlation between the galaxy infall time tinfall and age of the dynamically cold disk tcold for four different mass bins: 8 < log M*/M⊙ < 8.6 (light cyan), 8.6 < log M*/M⊙ < 9.2 (cyan), 9.2 < log M*/M⊙ < 10 (blue), and 10 < log M*/M⊙ < 12 (dark blue). Dashed lines mark the 1σ confidence limits of each correlation. Circles and crosses correspond to galaxies with extended and concentrated dynamically cold disks, respectively. Each galaxy is plotted by the age of its dynamically cold disk and infall time given by the median of the correlation corresponding to its stellar mass. The error bars corresponds to the 1σ uncertainty of the infall time inferred from the correlation. Right panel: distribution of the left-panel galaxies in the phase-space of the projected line-of-sight velocity of the galaxy normalized by line-of-sight velocity dispersion of the cluster VLOS/σLOS versus the projected clustercentric radius of the galaxy normalized by the cluster virial radius Rproj/Rvir. Each galaxy is color-coded according to its infall time except for the 4 LTGs without age and metallicity information shown in gray. The symbols are the same as in the left panel and the LTGs are marked by gray symbols. The boundaries of the regions A, B, C, D, and E are defined as in Rhee et al. (2017).

For galaxies with the cold-disk age obtained from the population-orbit superposition model, we can thus use this correlation to infer their infall time into the Fornax cluster. As shown in Fig. 5, we obtained the infall time of each galaxy by means of the correlation for the mass bin including its stellar mass. The values inferred from the median and 1σ scatter of the correlation are taken as the median and 1σ uncertainty of the galaxy infall time t_infall. They are listed in the last column of Table 2.

5.2. Comparison with galaxy locations in the projected phase-space

We cross-checked the infall time inferred from the cold-disk age with the galaxy location in the phase-space V_LOS/σ_LOS vs. R_proj/R_vir, as shown in the right panel of Fig. 5. Here, V_LOS/σ_LOS is the projected line-of-sight velocity of the galaxy normalized by line-of-sight velocity dispersion of the cluster and R_proj/R_vir is the projected clustercentric radius of the galaxy normalized by the cluster virial radius. The galaxies are color-coded by the infall time inferred from the cold-disk age in the left panel of Fig. 5.

Following Iodice et al. (2019b), we define ancient infallers as galaxies that fell into the cluster at 8–12 Gyr ago, intermediate infallers at 4–8 Gyr ago, and recent infallers at 0–4 Gyr ago. We divided the phase-space into five regions like in Iodice et al. (2019b), with the fraction of ancient infallers decreaseing from regions E and D to B+C (based on Rhee et al. 2017). About 50% of the galaxies in region E are assumed to be ancient infallers in our sample, while six out of seven galaxies in region E are ancient infallers and one is intermediate infaller. Galaxies in region D are supposed to be dominated by recent and intermediate infallers, while we have two ancient infallers and two recent infallers. Galaxies in region B+C are suggested to be dominated by recent infallers, while we have two ancient, one intermediate, and two recent infallers. The fractions of ancient, intermediate and ancient infallers in different regions we obtained and in comparison to that predicted from Rhee et al. (2017) are shown in Table 3.

Although we have a low amount of statistics at hand, we remain consistent with Rhee et al. (2017) in that the fraction of ancient infallers is highest in region E and then lower in D and B+C. However, the fraction of ancient infallers in our sample is generally higher, as we only included ETGs – and not the whole sample of galaxies in the Fornax cluster. The infall time inferred from cold-disk age seems reasonable for these Fornax cluster galaxies. Importantly, the correlation t_infall versus t_cold goes beyond a mere statistical estimation, providing a proxy of the cluster infall time per individual galaxy.

6. Dependence of cold disk assembly on galaxy mass and cluster environment

In this section, we show the dependence of the cold disk properties, including the luminosity fraction, age, and metallicity on galaxy stellar mass and cluster environment.

6.1. Cold-disk luminosity fraction

We studied the dependence of cold-disk fraction on galaxy stellar mass and infall time into the cluster. The cold-disk fraction is calculated as the fraction of cold orbits (λ_z ≥ 0.8) within 1R_e and 2R_e, respectively, of each galaxy.

We show the cold-disk fraction within 1 R_e in Fig. 6. In the left panel, we show the cold-disk fraction as a function of stellar mass. The 16 galaxies with colored models are color-coded by infall time inferred from their cold-disk age, while the four LTGs without age information in the models are considered as recent infallers, as they are dominated by younger stars at all radii.

Fig. 6. Dependence of cold-disk fraction (within 1Re) on galaxy stellar mass and infall time. Left panel: cold-disk fraction of the sample galaxies as a function of stellar mass, color-coded by infall time into the Fornax Cluster to highlight the ancient (red), intermediate (white), and recent (blue) infallers. Gray symbols correspond to the four LTGs without colored model which are considered as recent infallers. The gray curve represents the cold-disk fraction with 1Re of CALIFA galaxies (Zhu et al. 2018b). Circles and crosses correspond to galaxies with extended and concentrated dynamically cold disks, respectively. Right panel: cold-disk fraction within 1Re as a function of infall time color-coded by stellar mass, divided into three mass bins: log M*/M⊙ < 9.2 (light blue), 9.2 < log M*/M⊙ < 10 (blue), and log M*/M⊙ > 10 (dark blue). The four LTGs without colored model are plotted with an infall time upper limit of 4 Gyr ago.

Generally, the recent and intermediate infallers have higher cold-disk fractions, the galaxies with t_infall < 9 Gyr have cold-disk fractions, which vary as a function of stellar mass, comparable to the cold-disk fractions of field galaxies in the CALIFA survey (Zhu et al. 2018b). The most ancient infallers (red and orange points with t_infall > 9 Gyr) have low cold-disk fractions of f_cold(r < R_e)≲0.1 at all mass regions.

In the right panel, we show cold-disk fraction as a function of galaxy infall time by dividing the galaxies into three mass bins. For galaxies with similar stellar mass, the cold-disk fraction decreases from recent to ancient infallers. Although there is some scatter, especially in the case of FCC 153, with t_infall ∼ 9 Gyr but a high cold-disk fraction of f_cold(< 1R_e) > 0.3.

Similar figures but for cold-disk fraction within 2 R_e are shown in Fig. 7. The cold-disk fraction within 2 R_e is about 0.1 higher than that within 1 R_e for the galaxies with extended disks. The cold-disk fractions of f_cold(< 1R_e)∼0.3 and f_cold(< 2R_e)∼0.4 at M_∗ ∼ 10¹⁰ M_⊙ for the recent infallers are a factor of ∼3× or ∼4× higher than the ancient infallers mostly with both f_cold(< 1R_e) < 0.1 and f_cold(< 2R_e) < 0.1.

Fig. 7. Dependence of cold-disk fraction (within 2Re) on galaxy stellar mass and infall time. Details are similar to Fig. 6, but for cold-disk fractions within 2Re, which are about 0.1 higher than that within 1Re for galaxies with extended disks.

6.2. Cold-disk age and age gradient

We further studied the age and age gradient of cold disks in these galaxies. We calculated the age gradients of cold disks ∇t_cold for the nine galaxies with extended disks and colored models. As we can see from Fig. 4, the age profiles of disks are diverse in the inner regions and they become flat at the outer regions for most galaxies. We calculated the age gradients with profiles at r < R_e consistently for all galaxies, except for FCC 170 and FCC 255, the age gradient is calculated within the whole data coverage as their disks have very low surface-brightness at r < R_e.

In the left panel of Fig. 8, we show the mean age of the cold disk vs. stellar mass for the nine ETGs with extended disks and with colored model. The disks are generally older in the more massive galaxies. In the right panel of Fig. 8, we show the age gradient of the cold disk versus stellar mass colored by the galaxy infall time into the cluster. We get a positive or zero age gradient of the cold disk for most galaxies. The least-massive galaxies show significantly positive ∇t_cold for r < R_e and the age profiles become shallower at the outer regions. Three of the most massive galaxies FCC 083, FCC 147, and FCC 167 still have negative ∇t_cold at r < R_e, while their age profiles become flat or even turn to be positive (for FCC 167) at the outer regions.

Fig. 8. Mean age (left panel) and age gradient (right panel) of the extended cold disks in the sample galaxies as a function of stellar mass and infall time into the Fornax Cluster. Dashed line in the right panel marks the zero age gradient.

6.3. Cold-disk metallicity and metallicity gradient

In the left panel of Fig. 9, we show the mean metallicity of the cold disks versus stellar mass M_* colored by the galaxy infall time into the cluster. We only show the nine galaxies with extended cold disks and colored models. The cold disks appear to be similarly metal rich from low mass to high mass galaxies. The average metallicity of the cold disks is ⟨Z_cold⟩ = 0.91 ± 0.6 Z_⊙, with no significant correlation with galaxy stellar mass or infall time into the cluster. Comparing to the general metallicity-mass relation of galaxies (Gallazzi et al. 2005), the cold disks in galaxies with M_* ≲ 10¹⁰ M_⊙ in our sample are systematically more metal-rich than the general population of galaxies, while in more massive galaxies, they are more similar in composition.

Fig. 9. Mean metallicity (left panel) and metallicity gradient (right panel) of the extended cold disks in the sample galaxies as a function of stellar mass and infall time into the Fornax Cluster. The black curve and gray shadow represent the metallicity–stellar mass relation and the 1σ uncertainty for the galaxies in Gallazzi et al. (2005). Symbols are same as Fig. 8.

In the right panel of Fig. 9, we show the metallicity gradient of the cold disks vs. stellar mass. Metallicity gradients are calculated in a way similar to age gradients and only for the nine galaxies with extended disks and colored models. Almost all these galaxies show negative metallicity gradients for the cold disks, which are more metal-rich in the inner regions than in the outer one. There is no significant correlation seen between the metallicity gradient of cold disk and stellar mass – other than the change of metallicity gradients of the whole galaxy from positive to negative for low to high mass galaxies (Goddard et al. 2017; Zhuang et al. 2019), however, with low statistics of only nine galaxies included.

6.4. Cold versus hot galaxy components

We compare the age and metallicity of the dynamically cold disk and dynamically hot non-disk component in Fig. 10. The stars in the cold disks are as old as those of the hot components for the ancient infallers, whereas they are significantly younger than those of the hot components in recent infallers. The cold disks are in general more metal-rich that the hot components, which could explain the difference between our results on cold disk and the general mass-metallicity relation (Gallazzi et al. 2005), as shown in Fig. 9.

Fig. 10. Comparison of the mean ages (left panel) and mean metallicities (right panel) of the dynamically cold disks and dynamically hot non-disk components. Circles and crosses correspond to galaxies with extended and concentrated cold disks, respectively. The mean ages and metallicities are calculated for r < Re (filled circles). Black dashed lines represent y = x.

We compared the age and metallicity gradients of the two components in Fig. 11. For the galaxies with negative age gradients in cold disks, they also have negative age gradients in the hot component, while for the galaxies with positive age gradient in cold disks, the age gradients in their hot component could be very different; whereas both the cold disk and hot component have negative metallicity gradients in most galaxies and are generally consistent with each other.

Fig. 11. Comparison of the age gradient (left panel) and metallicity gradient (right panel) of the dynamically cold disks and dynamically hot non-disk components, for the galaxies with extended disks. Black dashed lines represent y = x. Blue dashed lines represent zero gradients.

The age and metallicity gradients of the cold disks are similar to that of the whole galaxy, as shown in Fig. 12. The signs of the gradients in cold disks are the same as what is seen for the whole galaxy in most cases. In this sense, the gradients directly measured for the whole galaxy could statistically reflect that in the cold disks, for galaxies with extended cold disks.

Fig. 12. Comparison of the age gradient (left panel) and metallicity gradient (right panel) of the dynamically cold disks and the whole galaxy, for the galaxies with extended disks. Symbols are the same as in Fig. 11.

7. Discussion

We find that the most ancient infallers in the Fornax cluster have low cold-disk fractions of f_cold ≲ 0.1, which is lower by about a factor of 4 than those of the recent and intermediate infallers in the Fornax cluster and of the field galaxies (Zhu et al. 2018b), with similar stellar mass. Although there is some scatter, especially with regard to FCC 153, which is potentially an ancient infaller, but with a high cold-disk fraction.

Our results indicate that a galaxy’s infall into a cluster has a strong effect on the assembly of the cold disk, which might be caused by the cut-off of the gas accretion from the cluster environment or removal of the gas by ram pressure and tidal stripping. Both processes explain the above results: galaxies that have earlier fallen into the cluster have their gas earlier removed or cut off, thus their disks are older and contribute a smaller fraction to total luminosity; galaxies that have recently fallen into the cluster have more time to grow their disk and the disk stars could be younger. However, it should be noted that we have the luminosity and not the mass fraction of the disk. The higher luminosity fraction of the cold disks in recent infallers could be partially caused by the lower stellar mass-to-light ratio in younger stellar disks.

It has been argued that ancient infallers could have a distinct evolutionary history before entering the cluster environment, which determined their morphology, kinematics, and gas content “at infall” (e.g., Han et al. 2018; Su et al. 2021). This could be part of the reason for the difference we see between the ancient and recent infallers at present time. But this preprocessing cannot explain all the differences. Numerical simulations show that dynamically cold disks are often partially removed or heated by tidal disruption in cluster (Joshi et al. 2020). With the control of galaxy properties at the time of infall, the galaxies that fall into a cluster result in significantly smaller disk fraction at redshift zero than those not. In this sense, the ancient infallers now in the cluster center would suffer most from the tidal heating or disruption, thus their cold-disk fraction could be further reduced.

We find positive or flat age gradients in the extended disks, which is opposite to what expected from a normal inside-out disk growth. The turnover of disk-age gradients could be the result of the continuous star formation in the galaxy inner regions, while the star formation in the outer regions is quickly reduced after falling into the cluster. This is generally consistent with the scenario found in the EAGLE simulations (Pfeffer et al. 2022). In five of our nine galaxies with extended disks, the stars in the inner disk region are about 2–5 Gyr younger than those in the outer disk region, which indicates that the star formation in galaxy inner regions could last for a few Gyrs more after the galaxy falls into the cluster. None of the 15 ETGs of our sample exhibit any detectable ongoing star formation: stars in the inner disks are younger, but their star formation stopped a few Gyrs ago. On the other hand, the inner regions are usually blended with an old bulge component, the age gradient of the whole galaxy is thus not so obviously negative as the disk only. Lasting star formation in galaxy inner regions is usually difficult to be directly observed in galaxies at present time, while hopefully might be directly observable at high redshift (Maier et al. 2019a).

Strong positive age gradients have been reported for the SAURON ETGs in the mass range of 10¹⁰⁻¹² M_⊙ using the method based on line indices (Kuntschner et al. 2010). However, in larger sample of MaNGA galaxies, the age gradients of ETGs are found to be generally shallow (Goddard et al. 2017). Positive age gradients are mostly found in low mass ETGs with M_* ≲ 10^10.5 (Li et al. 2018) by using the technique of full spectra fitting to derive the stellar populations. As shown in Li et al. (2018), the strong positive age gradients reported in Kuntschner et al. (2010) could disappear when uniformly using the same technique. As we are also using the full spectra fitting, a comparison with these recent results from MaNGA galaxies might offer a more fair outcome. We note that the age gradients found in MaNGA galaxies are different in Goddard et al. (2017) and Li et al. (2018), while the samples these authors studied are also different.

We found positive age gradients in the cold disks of least-massive galaxies, namely, M_* ≲ 10^10.5, whereas the disks that have negative age gradients in the inner regions and tend to be flat in the outer regions of most massive galaxies. This is consistent with the positive age gradient in low-mass MaNGA ETGs (Li et al. 2018), considering age gradients of the cold disks are statistically consistent with that of the whole galaxies. The difference of age gradients in low mass and high mass ETGs could be explained by the difference of gas content in low and high mass galaxies at infall. The low mass galaxies could have large gas fractions when they fall into the cluster. The gas in the outer regions is stripped, but the gas kept in inner regions form a significant amount of new stars that finally turn-over the age gradient. On the other hand, massive galaxies could have already formed most of their stars before entering the cluster. They have small gas fractions when falling into the cluster, thus the age gradient is less likely to be turned over by the formation of new stars after entering the cluster. Although with limited statistics, there is no galaxy with significant negative age gradient in our sample with M_* ≲ 10^10.5. Approximately flat gradients for such galaxies are found in the MaNGA sample studied by Goddard et al. (2017), while the MaNGA sample studied in Li et al. (2018) includes galaxies in different environment, and shows large scatter of age gradients of ETGs in the same mass range, which may indicate the effects of different environments. Although other physical processes, such as mergers, could also cause new star formation in the center of these relatively low mass ETGs in the field (Yozin & Bekki 2012; Lagos et al. 2022).

Similar ages of the disk and hot components in ancient infallers indicate that stars also formed on dynamically warmer orbits after falling. Star formation onto dynamically warmer orbits is possible if the gas is compacted into the inner regions of the galaxy, resulting in a high gas density in the galaxy center after falling into the cluster. As shown in Fig. 4, we also find younger stars in the very center of the hot components for a few galaxies, which is consistent with the aforementioned scenario.

8. Conclusions

We built orbit-superposition models for 20 galaxies (15 ETGs + 5 LTGs) in the Fornax cluster observed with MUSE/VLT in the context of the Fornax3D survey. We built a population-orbit superposition model for 16 sample galaxies by coloring the orbits with age and metallicity. By simultaneously fitting the surface-brightness, stellar kinematic, age, and metallicity maps, we obtain the internal stellar orbit distribution, as well as age and metallicity distribution of each galaxy. Based on the best-fit population-orbit superposition model, we decompose each galaxy into a dynamically cold disk (λ_z ≥ 0.8) and a dynamically hot non-disk component (λ_z < 0.8). Then, we obtain the surface-brightness, age, and metallicity radial profiles for each component by projecting them face-on. We analyze the dependence of the cold disk properties on galaxy stellar mass and infall time into the Fornax cluster. Our main results are as follows:

1. We estimate the galaxy infall time into the cluster based on a tight correlation with cold-disk age calibrated with the cosmological simulation TNG 50. The infall time inferred in this way is statistically consistent with the galaxy’s location in the cluster’s observational phase-space diagram.

2. The ancient infallers have significantly lower luminosity fractions of the cold disk component, regardless of stellar mass. The recent and intermediate infallers have cold-disk fraction increasing as a function of stellar mass. They reach f_cold(< R_e)∼0.3 and f_cold(< 2R_e)∼0.4 at M_∗ ∼ 10¹⁰ M_⊙, consistent with the CALIFA galaxies in the field. Most ancient infallers have f_cold ≲ 0.1, which is a factor of ∼4 lower than that of recent infallers of the same stellar mass.

3. Nine of the 16 galaxies with population-orbit superposition have extended cold disks. The five least-massive galaxies have positive or flat age gradients in their cold disks. Indeed, the stars in the inner disk regions are about 2–5 Gyr younger than those in the outer disk regions. In contrast, the four more massive galaxies, which might have a lower gas content at infall, have negative age gradients in the inner disk regions and flat age profiles in the outer disk regions.

Our results indicate that the star formation rate in disks can be significantly reduced after the galaxies fall into a cluster in general. While the star formation in the outer regions of a galaxy could stop shortly after it falls into the cluster, star formation in the inner regions should last for longer time causing the stars in the inner disk regions to be 2–5 Gyr younger than those in the outer disk regions. This is generally consistent with the scenario found for the Virgo cluster (Crowl & Kenney 2008), while we have been able for the first time to quantify the age difference between the inner and outer regions of cold disks.

By taking advantage of the high-quality MUSE datacubes and structure decomposition based on the population-orbit superposition models, we have been able to isolate and accurately study the dynamically cold disk component. This allows for an in-depth and direct comparison of the galactic structures with galaxies formed in cosmological simulations, which will help us to further understand the physical processes driving galaxy evolution.

Acknowledgments

The models and corresponding results presented in this paper would not have been possible without the VLT-MUSE data granted through the award of DDT telescope time (ESO programme 296.B-5054(A)). LZ acknowledges the support from the National Key R&D Program of China under grant No. 2018YFA0404501, National Natural Science Foundation of China under grant No. Y945271001, and CAS Project for Young Scientists in Basic Research under grant No. YSBR-062. GvdV acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 724857 (Consolidator Grant ArcheoDyn). LC acknowledges financial support from Comunidad de Madrid under Atracción de Talento grant 2018-T2/TIC-11612 and the Spanish Ministerio de Ciencia, Innovación y Universidades through grant PGC2018- 093499-B-I00. JF-B, IMN and FP acknowledge support through the RAVET project by the grant PID2019-107427GB-C32 from the Spanish Ministry of Science, Innovation and Universities (MCIU), and through the IAC project TRACES which is partially supported through the state budget and the regional budget of the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community. The F3D data is based on observations collected at the European Southern Observatory under ESO programme 296.B-5054(A), and available in the ESO Science Archive Facility. EMC is supported by MIUR grant PRIN 2017 20173ML3WW-001 and Padua University grants DOR2019-2021.

References

Appendix A: Mass density MGE of all galaxies

Appendix B: Best-fit models of all galaxies

Fig. B.1. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 083. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.2. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 119. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.3. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 143. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.4. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 147. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.5. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 148. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.6. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 153. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.7. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 161. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.8. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 167. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.9. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 170. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.10. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 179. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively. We find a possible weak bar structure at |X| ≲20 arcsec and |Y| ≲10 arcsec.

Fig. B.11. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 182. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.12. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 249. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.13. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 255. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.14. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 263. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.

Fig. B.15. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 276. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.16. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 290. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.

Fig. B.17. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 301. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.

Fig. B.18. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 308. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.

Fig. B.19. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 312. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.

Appendix C: Grid parameters of FCC 177

Fig. C.1. Optimization over the 5D parameter-space for FCC 177. Each points indicates an exploration of the parameter-space with the smaller orbit sampling. Color represents the χ2 of each parameter set. The best-fit parameter is indicated by the yellow most point with a cross symbol.

All Tables

All Figures

	Fig. 1. Best-fit population-orbit superposition model of FCC 177. From top to bottom: maps of the data, model, and residuals (model−data). From left to right: maps of the surface-brightness, SB, mean velocity, V, velocity dispersion, σ, Gauss-Hermite coefficients h3 and h4, light-weighted age and metallicity of the stars. Similar plots for the other sample galaxies are shown in Figs. B.1–B.19.
In the text

Fig. 2. Orbital decomposition of FCC 177. Probability density distribution p(r, λz) in the left panel, age distribution p(r, t) in the central panel, and metallicity distribution p(r, Z) in the right panel of the stellar orbits in the phase space of time-averaged radius, r, vs. circularity, λz. The probability densities are normalized to unity within the data coverage. All the distributions are averages of multiple best-fit models that fall within the 1σ confidence level of the model hyper-parameters. The dashed line marks our orbit-based division into two components: a dynamically cold disk component (λz ≥ 0.8) and a dynamically hot non-disk component (λz < 0.8). The shadow regions are beyond the data coverage of age and metallicity maps. Similar plots for the other galaxies are shown in Figs. B.1–B.19.

In the text

Fig. 3. Surface-brightness, age, and metallicity radial profiles of the whole galaxy, the dynamically cold disk, and the dynamically hot non-disk component for FCC 177. We show the profiles of the whole galaxy (left panels) with the blue solid curve, and the profiles of dynamically cold disk and dynamically hot non-disk component (right panels) with the blue solid curve and black dashed curve. The shadowed areas indicate the scatter of the profiles of models that fall within the 1σ confidence level. The age and metallicity profiles of the dynamically cold disks are considered reliable and shown in the regions where the dynamically cold disk contributes at least 10% of the total surface brightness.

In the text

Fig. 4. Surface-brightness, age, and metallicity radial profiles of the galaxy, dynamically cold disk, and dynamically hot non-disk component for the sample galaxies. We divided the galaxies into two groups: 11 galaxies with a dynamically cold disk extended out to large radii (left panels) and 9 galaxies with a dynamically cold disk concentrated in the inner regions (right panels). The radial profiles are color-coded by the galaxy stellar mass. The shadowed areas indicates the scatter of the profiles of models that fall within the 1σ confidence level of the model hyper-parameters. Dashed lines refer to the four LTGs without orbits tagged with age and metallicity. The age and metallicity radial profiles of the dynamically cold disks are considered reliable and shown in the regions where the dynamically cold disk contributes at least 10% of galaxy surface-brightness. Therefore, these profiles are available only for a limited radial range for the galaxies with a concentrated cold disk.

In the text

Fig. 5. Infall time into the Fornax cluster of the sample galaxies for which we obtained the age of the dynamically cold disk. Left panel: correlation between the galaxy infall time tinfall and age of the dynamically cold disk tcold for four different mass bins: 8 < log M*/M⊙ < 8.6 (light cyan), 8.6 < log M*/M⊙ < 9.2 (cyan), 9.2 < log M*/M⊙ < 10 (blue), and 10 < log M*/M⊙ < 12 (dark blue). Dashed lines mark the 1σ confidence limits of each correlation. Circles and crosses correspond to galaxies with extended and concentrated dynamically cold disks, respectively. Each galaxy is plotted by the age of its dynamically cold disk and infall time given by the median of the correlation corresponding to its stellar mass. The error bars corresponds to the 1σ uncertainty of the infall time inferred from the correlation. Right panel: distribution of the left-panel galaxies in the phase-space of the projected line-of-sight velocity of the galaxy normalized by line-of-sight velocity dispersion of the cluster VLOS/σLOS versus the projected clustercentric radius of the galaxy normalized by the cluster virial radius Rproj/Rvir. Each galaxy is color-coded according to its infall time except for the 4 LTGs without age and metallicity information shown in gray. The symbols are the same as in the left panel and the LTGs are marked by gray symbols. The boundaries of the regions A, B, C, D, and E are defined as in Rhee et al. (2017).

In the text

Fig. 6. Dependence of cold-disk fraction (within 1Re) on galaxy stellar mass and infall time. Left panel: cold-disk fraction of the sample galaxies as a function of stellar mass, color-coded by infall time into the Fornax Cluster to highlight the ancient (red), intermediate (white), and recent (blue) infallers. Gray symbols correspond to the four LTGs without colored model which are considered as recent infallers. The gray curve represents the cold-disk fraction with 1Re of CALIFA galaxies (Zhu et al. 2018b). Circles and crosses correspond to galaxies with extended and concentrated dynamically cold disks, respectively. Right panel: cold-disk fraction within 1Re as a function of infall time color-coded by stellar mass, divided into three mass bins: log M*/M⊙ < 9.2 (light blue), 9.2 < log M*/M⊙ < 10 (blue), and log M*/M⊙ > 10 (dark blue). The four LTGs without colored model are plotted with an infall time upper limit of 4 Gyr ago.

In the text

	Fig. 7. Dependence of cold-disk fraction (within 2Re) on galaxy stellar mass and infall time. Details are similar to Fig. 6, but for cold-disk fractions within 2Re, which are about 0.1 higher than that within 1Re for galaxies with extended disks.
In the text

	Fig. 8. Mean age (left panel) and age gradient (right panel) of the extended cold disks in the sample galaxies as a function of stellar mass and infall time into the Fornax Cluster. Dashed line in the right panel marks the zero age gradient.
In the text

	Fig. 9. Mean metallicity (left panel) and metallicity gradient (right panel) of the extended cold disks in the sample galaxies as a function of stellar mass and infall time into the Fornax Cluster. The black curve and gray shadow represent the metallicity–stellar mass relation and the 1σ uncertainty for the galaxies in Gallazzi et al. (2005). Symbols are same as Fig. 8.
In the text

	Fig. 10. Comparison of the mean ages (left panel) and mean metallicities (right panel) of the dynamically cold disks and dynamically hot non-disk components. Circles and crosses correspond to galaxies with extended and concentrated cold disks, respectively. The mean ages and metallicities are calculated for r < Re (filled circles). Black dashed lines represent y = x.
In the text

	Fig. 11. Comparison of the age gradient (left panel) and metallicity gradient (right panel) of the dynamically cold disks and dynamically hot non-disk components, for the galaxies with extended disks. Black dashed lines represent y = x. Blue dashed lines represent zero gradients.
In the text

	Fig. 12. Comparison of the age gradient (left panel) and metallicity gradient (right panel) of the dynamically cold disks and the whole galaxy, for the galaxies with extended disks. Symbols are the same as in Fig. 11.
In the text

	Fig. B.1. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 083. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.2. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 119. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.3. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 143. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.4. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 147. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.5. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 148. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.6. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 153. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.7. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 161. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.8. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 167. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.9. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 170. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.10. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 179. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively. We find a possible weak bar structure at |X| ≲20 arcsec and |Y| ≲10 arcsec.
In the text

	Fig. B.11. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 182. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.12. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 249. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.13. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 255. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.14. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 263. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.
In the text

	Fig. B.15. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 276. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.16. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 290. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.
In the text

	Fig. B.17. Best-fit population-orbit superposition model (top panels) and orbital decomposition (bottom panels) of FCC 301. The content of the top and bottom panels is the same as in Figs. 1 and 2, respectively.
In the text

	Fig. B.18. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 308. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.
In the text

	Fig. B.19. Best-fit orbit superposition model (left panels) and orbital decomposition (right panels) of FCC 312. Details are similar to Figs. 1 and 2 as we show for FCC 177, but we do not show the age and metallicity.
In the text

	Fig. C.1. Optimization over the 5D parameter-space for FCC 177. Each points indicates an exploration of the parameter-space with the smaller orbit sampling. Color represents the χ2 of each parameter set. The best-fit parameter is indicated by the yellow most point with a cross symbol.
In the text