© L. Zhu et al. 2022

1. Introduction

Driven by gravity, galaxies are expected to continuously grow through the merging of smaller systems. Deriving their past merger history is challenging, as the accreted stars disperse quickly. Uncovering the properties of past merger events has so far only been possibly in the nearest few galaxies, where individual stars can be resolved in observations.

In our own galaxy, the Milky Way, the motion and chemistry of individual stars are detected. Numerous streams and over-densities (e.g., Majewski et al. 2003; Helmi 2020; Ibata et al. 2021) identified in the outer halo indicate the many minor mergers still ongoing in the Milky Way. Debris from massive ancient mergers is spatially well mixed, but fortunately still distinguishable in kinematics. Recently, stars on highly radial orbits in the vicinity of the Sun have been associated with the Galaxy’s inner stellar halo – for example, the so-called Gaia Enceladus or Gaia Sausage, which is thought to have originated from a single massive merger with stellar mass of ∼3−6 × 10⁸ M_⊙ about 10 Gyr ago (Helmi et al. 2018; Belokurov et al. 2018).

In our neighboring galaxy, the Andromeda giant spiral galaxy (M 31), an inner stellar halo with a giant stream has been identified (Ibata et al. 2014). It could also have been induced by a single massive merger with M 32p that occurred about ∼2 Gyr ago (D’Souza & Bell 2018a; Hammer et al. 2018) with a stellar mass of ∼2 × 10¹⁰ M_⊙. The satellite mass is inferred from the metallicity of M 32, which is supposed to be the surviving core of the progenitor satellite galaxy. The merger time was inferred from the stellar age distribution of stars in the giant streams and inner stellar halo. In the nearest giant elliptical galaxy, NGC 5128 at distance of 3.8 Mpc, individual stars in the outer halo are resolved by deep Hubble Space Telescope photometry (Rejkuba et al. 2011), and analysis similar to that of M 31 indicates a past massive merger of also ∼2 × 10¹⁰ M_⊙, (D’Souza & Bell 2018b).

In the above cases, the quantification of the satellite mass relies on the mass-metallicity relation, while the constraint on merger time relies on star formation histories obtained from individual stars. Similar studies are thus limited to the nearest few galaxies until 30 m telescopes become available and allow us to obtain resolved stellar populations for more distant galaxies. The past merger histories for a few galaxies at larger distances have been partially uncovered using different techniques. Individual globular clusters in the halos are detected in nearby galaxies and are used as another tracer of the galaxy assembly history (Forbes et al. 2016; Beasley et al. 2018; Kruijssen et al. 2019). Deep images broadly characterize the light distribution of stellar halos in nearby galaxies, and they have been used to find various substructures, such as streams, in the halo (Merritt et al. 2016; Harmsen et al. 2017).

Over the last two decades, integral-field-unit (IFU) instruments have spectroscopically mapped thousands of galaxies across a wide range of masses and Hubble types (de Zeeuw et al. 2002; Cappellari et al. 2011; Sánchez et al. 2012; Bryant et al. 2015; Bundy et al. 2015). In principle, information regarding stellar motions and chemical distributions of an external galaxy is included in these IFU data, though all mixed along the line-of-sight. Recent studies have attempted to constrain the global ex situ fractions (Davison et al. 2021a,b; Boecker et al. 2020) or the mass of satellite mergers (Pinna et al. 2019a,b; Martig et al. 2021) based on age and metallicity distributions of stars in galaxies’ inner regions obtained from IFU data. These methods can identify minor mergers as their accreted stars have lower metallicity than in situ stars, but they become insensitive to major mergers with overlapping metallicities. Also, the merger times remain largely unconstrained with such approaches.

From a numerical perspective, recent and current cosmological hydrodynamical simulation projects, such as IllustrisTNG¹ and EAGLE², have been able to reproduce large numbers of galaxies across the morphological spectrum with well-resolved structures (Pillepich et al. 2019; Pulsoni et al. 2020; Rodriguez-Gomez et al. 2019; Du et al. 2019; Correa et al. 2017). Structures such as Gaia-Enceladus in the Galactic inner stellar halo also arise in Milky Way-like simulated galaxies (e.g., Grand et al. 2020). In Zhu et al. (2022), a hot inner stellar halo, consisting of orbits that are dominated by random instead of ordered motion, is defined universally for all types of galaxies across a wide mass range. The mass of such a hot inner stellar halo is found to be strongly correlated with the stellar mass of the most massive merger as well as with the total ex situ stellar mass. The correlations are consistent across the TNG50 (Pillepich et al. 2019; Nelson et al. 2019a), TNG100 (Pillepich et al. 2018; Nelson et al. 2018; Springel et al. 2018; Marinacci et al. 2018; Naiman et al. 2018), and EAGLE (Schaye et al. 2015; Crain et al. 2015) simulations for galaxies with stellar mass ≳4 × 10¹⁰ M_⊙, regardless of the numerical resolution and galaxy formation model. The hot inner stellar halo is thus a promising way to quantify the past merger mass for all types of massive galaxies.

In parallel, we have been able to accurately infer the internal stellar orbit distribution of external galaxies by applying Schwarzschild’s orbit-superposition dynamical method to the stellar kinematic maps, extracted from IFU data (van den Bosch et al. 2008; Zhu et al. 2018a). The method has been applied to hundreds of galaxies across the Hubble sequence (Zhu et al. 2018b,c; Jin et al. 2020; Santucci et al. 2022). Stellar orbits in the model have been further colored with age and metallicities in a recently developed population-orbit superposition method (Poci et al. 2019; Zhu et al. 2020). By fitting the stellar kinematic, age, and metallicity maps simultaneously, we have been able to obtain the internal stellar orbit distribution as well as the age and metallicity distributions (Poci et al. 2021). With this novel method, chemodynamical identification of galactic structures as employed in the Milky Way has become possible for external galaxies.

In this paper we show that the mass of the hot inner stellar halo is a quantity that can be robustly obtained for external galaxies by applying our population-orbit superposition method to Very Large Telescope (VLT)/MUSE IFU data of two Fornax cluster galaxies, NGC 1380 and NGC 1427, observed as part of the Fornax3D project (Sarzi et al. 2018). We subsequently contrast these robust mass measurements with the correlations obtained from simulations as presented in Zhu et al. (2022) to infer the masses of the mergers that induced the hot inner stellar halos in both galaxies. Finally, we use the modeled age distributions of the different dynamical structures to also constrain the times of the last major mergers in NGC 1380 and NGC 1427.

The paper is organized as follows. In Sect. 2 we describe the observational data of NGC 1380 and NGC 1427. In Sect. 3 we present the population-orbit superposition models of the two galaxies and measure the mass of their hot inner stellar halos. In Sect. 4 we present the main results of the paper and, after discussing effects due to the cluster environment in Sect. 5, we summarize in Sect. 6. Tests of our method and the evaluation of model uncertainties are outlined in Appendices A and B.

2. Deep imaging and spectroscopic data

2.1. NGC 1380 and NGC 1427 in the Fornax cluster

Astronomers have observed galaxies in the Fornax cluster in great detail in recent years. The Fornax Deep Survey (FDS) has observed all galaxies in the 9 square degrees around the core of the Fornax cluster with the VLT Survey Telescope down to 27 mag arcsec⁻² in the r band (Iodice et al. 2019a). Then, the Fornax3D project (Sarzi et al. 2018), using the MUSE IFU instrument on the VLT, observed all 23 early-type galaxies (ETGs) and 10 late-type galaxies within the virial radius of the cluster, down to 25 mag arcsec⁻² in the B band. NGC 1380 (FCC 167) and NGC 1427 (FCC 276) are two galaxies in the Fornax cluster observed as part of the FDS (Iodice et al. 2019a) and Fornax3D survey (Sarzi et al. 2018).

NGC 1380 is a lenticular galaxy located on the northwestern side of the cluster at ∼40 arcmin from the brightest central galaxy (BCG) NGC 1399, and at a distance of 21.2 Mpc from us (Blakeslee et al. 2009). NGC 1380 has a total r-band magnitude of m_r = 9.27 and a half-light radius of R_e = 56.2″ = 5.8 kpc. Deep imaging reveals an outer and fainter isophotes at R >  3 arcmin, which appear less flattened and twisted with respect to the inner and “diskier” regions (Iodice et al. 2019a). Through the Fornax3D project, NGC 1380 was observed with three MUSE pointings out to a galactocentric radius of 150″ along the major axis and up to 40″ along the minor axis (Fig. 1). It is one of the most-studied objects within the cluster. Some initial dynamical modeling results have been included in Sarzi et al. (2018).

Fig. 1. Population-orbit model of NGC 1380. Panel a: deep r-band image of NGC 1380 from the FDS with the VLT Survey Telescope (Iodice et al. 2019a). The three squares indicate the three pointings, each with a 1′ × 1′ square field-of-view, of spectroscopic observations with the VLT MUSE instrument. The black dashed ellipse corresponds to the isophote at μB = 25 mag arcsec−2. Panel b: data from the VLT/MUSE observations (first column), the best fitting population-orbit model (second column), and the residuals of the data minus model divided by the observational errors (third column). From top to bottom: surface brightness, mean velocity, velocity dispersion, GH coefficients h3 and h4, and the light-weighted age and metallicity maps of the stars in NGC 1380. We note that kinematic maps were derived from spectra binned with a different scheme than the age and metallicity maps.

NGC 1427 is the brightest elliptical galaxy located on the eastern side of the cluster at ∼50 arcmin from the BCG, and at a distance of 19.6 Mpc from us (Blakeslee et al. 2009). NGC 1427 has m_r = 10.15 mag and R_e = 44.7″ = 4.2 kpc. This galaxy becomes bluer with increasing radius and in the outskirts, and its outer and fainter isophotes appear asymmetric at R ≳ 3 arcmin (Iodice et al. 2019a). NGC 1427 was observed with two MUSE pointings out to a radius of 100″ along the major axis and up to 40″ along the minor axis (Fig. 2).

Fig. 2. Population-orbit model of NGC 1427. See the caption of Fig. 1 for a description.

Both NGC 1380 and NGC 1477 are located beyond the visible extended stellar halo of the BCG NGC 1399, which is mapped out to 33 arcmin (Iodice et al. 2016).

2.2. Stellar kinematic, age, and metallicity maps

The IFU spectra were spatially binned using a Voronoi tessellation (Cappellari & Copin 2003) to reach a minimal signal-to-noise S/N = 100 per bin, resulting in a total of 4300 bins for NGC 1380 and 3072 bins for NGC 1427. Following earlier analysis of NGC 1380 (Sarzi et al. 2018), the stellar kinematics were extracted by applying the penalized pixel fitting (pPXF; Cappellari & Emsellem 2004) to the spectra with wavelength range 4750−5500 Å. This yields 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₄.

To obtain high-quality maps of stellar age and metallicity, the spectra were spatially re-binned to reach S/N = 200 for NGC 1380, resulting in a total of 974 bins. For NGC 1427, we kept the 3072 bins with target S/N = 100 from above. Each binned spectrum was fitted with regularized pPXF (Cappellari 2017) using the MILES single stellar population models (Falcón-Barroso et al. 2011) based on BaSTI isochrones (Vazdekis et al. 2015), and we assume a constant Milky Way-like Kroupa (Kroupa 2002) stellar initial mass function. The models cover a wavelength range of 3540−7410 Å, and are sampled at a spectral resolution of 2.51 Å (Falcón-Barroso et al. 2011). They include 8 values of total metallicity [M/H] = [−0.96, +0.40] dex (with a resolution between 0.14 and 0.48 dex), 33 values of age between 1.0 and 14.0 Gyr (with a resolution between 0.25 and 0.5 Gyr), and two values of [Mg/Fe] = 0.0 (solar abundance) and 0.4 dex (super-solar). The templates with metallicity below −0.96 and age below 1.0 Gyr were excluded. We fitted the stellar kinematics and populations in the same run, using a multiplicative polynomial of eight order and a regularization parameter of 0.5. We generally followed the same approach as in earlier work (Pinna et al. 2019a) but renormalized the templates to get the light-weighted stellar age, t, and metallicity, Z/Z_⊙, in each Voronoi bin.

3. Population-orbital superposition method

We fit population-orbit superposition models to the luminosity distribution, stellar kinematic, age, and metallicity maps for each galaxy. The major steps of constructing the model include: (1) construction of the gravitational potential, (2) calculation of an orbit library in the gravitational potential, (3) fitting to the luminosity density and kinematic maps by weighting the orbits, and (4) fitting to the age and metallicity maps by coloring the orbits. Best fitting models are obtained by exploring the free parameters in the gravitational potential and minimizing the χ² difference between models and data. In the end, we obtain 3D models superposed by stellar orbits tagged with age and metallicity, which match all data well. We then decompose the galaxies based on the modeled stellar orbit distribution.

While our method was introduced and extensively verified in Zhu et al. (2020), we present the main steps below.

3.1. Gravitational potential

The gravitational potential is assumed to be generated by the stellar mass plus dark matter halo. Also, a central black hole of 10⁸ M_⊙, is added, although it remains unresolved by the kinematic data and does not affect our results. To obtain the stellar mass distribution, we first fit the r-band image from FDS by a multiple Gaussian expansion (MGE) model (Cappellari 2002). Then by adopting a set of viewing angles (ϑ, ψ, ϕ), we de-project toward a triaxial MGE model that represents the intrinsic stellar luminosity density. Finally, by multiplying with a constant stellar mass-to-light ratio M/L_r, we arrive at the intrinsic stellar mass distribution. The effects of the constant stellar mass-to-light ratio on our results are discussed in Appendix B.

The three viewing angles relate directly to three axis ratios (p, q, u) describing the intrinsic shape (van den Bosch et al. 2008). We leave the intrinsic intermediate-to-major and minor-to-major axis ratios (p, q) free to vary, but limit u in the range of 0.98−0.9999. Moderate triaxiality is allowed in the model. The dark halo is throughout assumed to be spherical with a NFW (Navarro et al. 1997) radial profile with two free parameters: the halo (virial) mass M₂₀₀, defined as the total mass enclosed within the virial radius r₂₀₀ within which the average density is 200 times the critical density, and the concentration c, defined as the ratio between dark matter virial radius and the dark matter scale radius. We let M₂₀₀ free and for each dark matter halo we fix c following the relation from Dutton & Macciò (2014). In total we thus have five free so-called hyperparameters: M/L_r, p, q, u, and M₂₀₀.

3.2. Fitting the data

3.2.1. Fitting to luminosity density and stellar kinematic maps

For each set of hyperparameters, we computed tens of thousands of orbits in the corresponding gravitational potential. The orbit sampling and calculation follow what we described in Zhu et al. (2020). The weights of the orbits in the resulting orbit library are then determined from fitting simultaneously the orbit-superposition to the intrinsic and projected luminosity density and stellar kinematic maps.

We adopted a parameter grid with intervals of 0.2, 0.02, 0.02, 0.01 and 0.2 for M/L_r, p, q, u, and log₁₀(M₂₀₀). Then we perform an iterative process searching for the best fitting models generally following the process described in Zhu et al. (2018a). As the initial searching steps were already small, we did not refine the search grid around the minimum. At the end, we got 294 models for NGC 1380 and 465 models for NGC 1427 with different collections of the hyperparameters. Once we explored the hyperparameter space, we selected the best fitting models with $\mathrm{\Delta }{\chi }^2\equiv {\chi }^2-min({\chi }^2)<\sqrt{2\times n_{\mathrm{GH}}\times N_{\mathrm{obs}}}$, where n_GH = 4 is the number of stellar kinematic moments and N_obs is the number of bins of each kinematic map. The criterion of $\sqrt{2\times n_{\mathrm{GH}}\times N_{\mathrm{obs}}}$ is roughly the χ² fluctuation of a single model when perturbing the kinematic data, it generally works well as the 1σ confidence level according to our test with mock data (Zhu et al. 2018a).

There are 55 models with different values of the hyperparameters within this 1σ confidence level for NGC 1380. We obtain an inclination angle of 77° ±2°, which is consistent with that estimated from the central dust ring (Viaene et al. 2019). Our dynamical model constrains the stellar mass of NGC 138 to be (1.33 ± 0.1) × 10¹¹ M_⊙, within 2 R_e, with M/L_r = 2.7 ± 0.3 M_⊙/L_⊙. The mass outside the kinematic data coverage is not well constrained. Therefore, we just take two times of the mass within R_e as NGC 1380’s total stellar mass, which gives M_* = (1.8 ± 0.2) × 10¹¹ M_⊙.

There are 150 models within this 1σ confidence level for NGC 1427. Our dynamical model gives a total stellar mass of (4.1 ± 0.3) × 10¹⁰ M_⊙, within 2 R_e, with M/L_r = 2.5 ± 0.2 M_⊙/L_⊙, and an inclination angle of 61° ±6°. Also in this case, the mass outside the kinematic data coverage is not well constrained, and so we take two times of the mass within R_e as NGC 1427’s total stellar mass, which gives M_* = (5.6 ± 0.6) × 10¹⁰ M_⊙.

3.2.2. Fitting to the age and metallicity maps

Next, for each of the models within the 1σ confidence level, we applied Voronoi binning in the phase space of r versus λ_z to group the orbits into ∼150 bundles with comparable weight. Then we tagged each bundle k with a single stellar age t_k and metallicity Z_k and compared the superposition of the orbit bundles with the light-weighted age and metallicity maps extracted from the VLT/MUSE observations.

We used Bayesian statistical analysis (Python package pymc3) to estimate the best fitting age and metallicity for each orbital bundle. To use the Bayesian theorem to compute the posterior of a model, we needed the prior and the data likelihood. We followed Zhu et al. (2020) in adopting bounded normal and log-normal distributions for the priors on t_k and Z_k. However, we adopted a Student’s t-distribution for the data likelihood (Salvatier et al. 2016), which will allow some outliers in the data and results in a robust fitting. For each parameter t_k (Z_k), we ran a chain with 2000 steps and took the last 500 steps to calculate the mean and 1σ values of t_k (Z_k) for each orbit bundle. Metallicity was fitted separately from age without implementing an age-metallicity correlation. The combination of the best fitting orbital weights with the best fitting age and metallicity of orbital bundles yields population-orbital models that fit in detail the observed stellar kinematic, age, and metallicity maps. The best fitting models of NGC 1380 and NGC 1427 are shown in Figs. 1 and 2, respectively.

3.3. Orbit-based decomposition of galaxies

Our method yields the weights of the different stellar orbits that contribute to the best fitting model as well as the average stellar age and metallicity of the stars represented by each orbit. We then characterize each stellar orbit with two main properties (Zhu et al. 2018a): the time-averaged radius, r, which represents the extent of the orbit, and the circularity, λ_z = J_z/J_max(E), which represents the angular momentum of the orbit around the minor z-axis normalized by the maximum of a circular orbit with the same binding energy E. Whereas |λ_z|∼1 represents dynamically cold orbits dominated by regular rotation, λ_z ∼ 0 represents dynamically hot orbits, which are dominated by radial random motions or long-axis tubes with regular rotation around the long axis. Negative values λ_z <  0 refer to orbits that counter-rotate with respect to the net (prograde) rotation. Figure 3 presents the resulting stellar-orbit distributions of NGC 1380 and NGC 1427, respectively, as probability-density distributions p(r, λ_z) in the phase space of r versus λ_z, as well as the age t(r, λ_z) and metallicity Z(r, λ_z) distributions of the stellar orbits in the same phase space.

Fig. 3. Orbital decompositions of NGC 1380 and NGC 1427. Panel a-1: probability density distribution of the stellar orbits p(r, λz) in the phase space of time-averaged radius r versus circularity λz in units of stellar mass per unit area in the phase space. We show the model within 140 arcsecs (or ∼14 kpc) for NGC 1380. The stellar masses are normalized such that they sum to unity within the data coverage. Panels b-1 and c-1: stellar age t(r, λz) and metallicity Z(t, λz) distribution of the orbits in the same phase space. All distributions p(r, λz), t(r, λz), and Z(t, λz) are averages of multiple best fitting models that fall within the 1σ confidence level of the model hyperparameters. The dashed lines indicate our orbit-based division into four components: disk, warm, bulge, and hot inner stellar halo. Panels a-2, b-2, c-2: same but for NGC 1427, which does not have an extended disk.

Considering the substructures in the stellar-orbit probability density and in the age and metallicity distributions in the phase space of r versus λ_z, we separated each galaxy into four stellar components: (1) disk (λ_z >  0.8), (2) bulge (λ_z <  0.8, r <  r_cut), (3) warm (0.5 <  λ_z <  0.8, r_cut <  r <  2 R_e), and (4) hot inner stellar halo (λ_z <  0.5, r_cut <  r <  2 R_e). The last is in fact a dynamically hot component, separated from the bulge by r_cut = 3.5 kpc, which is generally associated with a transition in density and metallicity distributions of the dynamically hot orbits.

In Fig. 4 we show the probability density and the average metallicity of the dynamical hot (λ_z <  0.5) and warm (0.5 <  λ_z <  0.8) orbits as a function of radius, for NGC 1380 (panels a-1 and b-1) and NGC 1427 (panels a-2 and b-2). For NGC 1380, both the probability density p_hot and metallicity Z_hot of the hot orbits first decreases sharply from an inner peak, and then flatten beyond a transition radius at ∼3.5 kpc. For NGC 1427, the transition is not so clear in the probability density, but there is a transition in the metallicity profile that flattens beyond ∼2 kpc. A transition generally exists in probability density profile of the hot orbits at ∼3.5 kpc for TNG50 galaxies, but especially those simulated galaxies with high ex situ fractions do not show a sharp transition similarly to NGC 1427 (Zhu et al. 2022). To be able to apply the resulting correlations obtained in (Zhu et al. 2022) for simulated galaxies, we adopt the same r_cut = 3.5 kpc to separate the bulge and hot inner stellar halo.

Fig. 4. Probability density and metallicity profiles. Panel a-1: probability density of dynamically hot orbits with λz <  0.5 (red) and warm orbits with 0.5 <  λz <  0.8 (orange) as a function of radius for NGC 1380. The red and orange thick and dashed curves are the mean and 1σ scatter of our best fitting models. The vertical line at r = 3.5 kpc marks the separation of the bulge from the hot inner stellar halo. Panel b-1: metallicity distribution of hot and warm orbits as a function of radius for NGC 1380. Panels a-2 and b-2: same but for NGC 1427.

The probability density p_warm and the average metallicity Z_warm of the dynamical warm orbits are also shown for comparison. The warm orbits are more metal-rich than hot orbits in NGC 1427, while the difference is small in NGC 1380. In both galaxies, the warm orbits show negative metallicity gradients from the inner to outer regions, which are different from the flattening of metallicity distribution of the hot orbits.

The four components are defined in such a way as to highlight distinct structural, kinematic, and stellar population distributions, thus likely implying distinct formation pathways. According to previous observational and numerical studies, it would appear that, for galaxies with a quiescent merger history, stars that formed in situ from the collapse of dense gas at high redshift constitute a compact spheroid (Du et al. 2021; Williams et al. 2013), which constitute a major part of our present-day “bulge” and part of the “warm component” (Yu et al. 2021). At lower redshift, in situ stars form out of a gaseous disk that subsequently leads to a dynamical cold stellar disk with strong rotation (Du et al. 2021; Yu et al. 2021). No substantial hot inner stellar halo is expected in galaxies with a quiet merger history (Zhu et al. 2022). On the other hand, according to the results of state-of-the-art cosmological galaxy simulations, the hot inner stellar halo is the product of massive mergers (Zhu et al. 2022), and stars in it can have three origins: (1) stars accreted from a massive satellite typically favor highly radial orbits with λ_z ∼ 0 (Boylan-Kolchin et al. 2008; Fattahi et al. 2019) and hence partially settle into the bulge and partially stay at larger radii, becoming part of the hot inner stellar halo; (2) at the same time, a merger event can destroy the previously existing disk or bulge of the main progenitor, inducing a fraction of the disk or bulge stars to be kicked to the phase-space regions of the hot inner stellar halo; and (3) a merger can last for a few billion years until coalescence and new stars formed during the merger event also partially contribute to the hot inner stellar halo.

As tested with mock data created from cosmological galaxy simulations, the luminosity fractions of such four components can be recovered in observed galaxies to good accuracy by our modeling (see Appendix A). In the cases of the Fornax cluster galaxies NGC 1380 and NGC 1427, the hot inner stellar halo contributes 27% and 34%, respectively, to the total r-band luminosity within 2 R_e. We estimated that the hot inner stellar halo mass is M_{*, halo(r <  2 Re)} = (3.6 ± 0.5) × 10¹⁰ M_⊙, and M_{*, halo(r <  2 Re)} = (1.4 ± 0.2) × 10¹⁰ M_⊙, for NGC 1380 and NGC 1427, respectively. The uncertainty includes the contribution of systematically errors evaluated based on rigorous tests with mock data and statistical errors of the model; the details are included in Appendix B.

4. Weighing and timing of past mergers

4.1. Stellar mass of the accreted satellite(s) by NGC 1380 and NGC 1427

Cosmological galaxy simulations show that the hot inner stellar halo is smoking-gun evidence of massive mergers that a galaxy has experienced in its past history. In a companion paper (Zhu et al. 2022), we show that cosmological galaxy simulations return strong correlations between the hot inner stellar halo mass M_{*, halo(r <  2 Re)} of galaxies and the stellar mass from the most massive galaxy they have ever merged with, M_*, Ex1. Combining all the stars from all accreted satellites and mergers, M_*, ExSitu, makes the correlations even stronger. These results hold for a definition of the hot inner stellar halo that is identical to the one put forward above and applied to the case of NGC 1380 and NGC 1427 via our population-orbit superposition models. The hot inner stellar halo mass M_{*, halo(r <  2 Re)} can thus be used to infer the mass of the most massive satellite (now destroyed) that a galaxy has ever merged with, as well as of the total accreted or ex situ stellar mass of an observed galaxy. This can be done thanks to the fact that the hot inner stellar halo mass can be well recovered by our population-orbit based model for real galaxies, even long after the stars from the merger events are phase-mixed (see Appendix A).

4.1.1. Stellar mass of the most massive satellite

We inferred the stellar mass of the most massive merger for NGC 1380 and NGC 1427 by using the scaling relations predicted by the cosmological simulation TNG50 (Pillepich et al. 2019; Nelson et al. 2019a) (see Fig. 5). For NGC 1380 with M_* = 1.8 × 10¹¹ M_⊙, and R_e = 5.8 kpc, we compared with the subsample of TNG50 galaxies with 8 × 10¹⁰ <  M_* <  4 × 10¹¹ M_⊙, and 4 <  R_e <  8 kpc. From these TNG50 galaxies, we obtained log₁₀(M_*, Ex1)=1.195 × log₁₀(M_{*, halo(r <  2 Re)})−2.04 with a Pearson correlation coefficient ℛ = 0.89, and a 1σ scatter against of the linear fit of σ = 0.23. With M_{*, halo(r <  2 Re)} = (3.6 ± 0.5) × 10¹⁰ M_⊙, for NGC 1380, we obtained M_*, Ex1 = 10^{10.57 ± 0.24} M_⊙, ($3.7_{-1.5}^{+2.7}\times {10}^{10}$ M_⊙), with a 3σ significant lower limit of 0.7 × 10¹⁰ M_⊙. The uncertainty of satellite mass is taken as $d=\sqrt{d_1^2+d_2^2}$, where d₁ is the uncertainty in the derivation of log₁₀(M_{*, halo(r <  2 Re)}), while d₂ is the uncertainty caused by the scatter of the relation (i.e., d₂ = σ).

Fig. 5. Stellar masses of the most massive satellites (now-destroyed) that merged with NGC 1380 and NGC 1427. Panel a: correlation from cosmological simulations between the mass of the hot inner stellar halo, M*, halo(r <  2 Re), and the stellar mass, M*, Ex1, of the most massive merged satellite. The black dots are a subsample of TNG50 galaxies with similar mass and size as NGC 1380, the thick black line is a linear fit to the points, and the thin black lines indicate the ±1σ scatter. The Pearson correlation coefficients, ℛ, between M*, halo(< 2 Re) and M*, Ex1 and 1σ scatter against the linear fit for the TNG50 subsample are indicated. The vertical red line with shadows marks the mass and corresponding ±1σ uncertainty of NGC 1380’s hot inner stellar halo. The horizontal dashed line marks the satellite(s) mass inferred from the correlations. The thick and thin curves in gray, blue, and green in the background are the running median and ±1σ scatter of the whole sample of TNG50, TNG100, and EAGLE galaxies. Panel b: same to panel a, but for the subsample of TNG50 galaxies with similar mass and size to NGC 1427. The vertical red line marks the mass of NGC 1427’s hot inner stellar halo, and the horizontal dashed line marks the mass of the satellite(s) inferred from the correlations.

In a similar fashion, for NGC 1427, with M_* = 5.6 × 10¹⁰ M_⊙, and R_e = 4.2 kpc, we compared with the subsample of TNG50 galaxies with 4 × 10¹⁰ <  M_* <  10¹¹ M_⊙, and 3 <  R_e <  7 kpc. From these TNG50 galaxies, we obtained log₁₀(M_*, Ex1)=1.21 × log₁₀(M_{*, halo(r <  2 Re)})−2.07 with the Pearson correlation coefficient ℛ = 0.77 and σ = 0.31. With M_{*, halo(r <  2 Re)} = (1.4 ± 0.2) × 10¹⁰ M_⊙, for NGC 1427, we inferred M_*, Ex1 = 10^{10.19 ± 0.32} M_⊙ ($1.5_{-0.7}^{+1.6}\times {10}^{10}$ M_⊙), with a 3σ significant lower limit of 0.2 × 10¹⁰ M_⊙.

4.1.2. Total accreted stellar masses

The correlations with M_{*, halo(r <  2 Re)} are stronger when combing all the ever accreted stellar mass M_*, ExSitu, as shown in Fig. 6. The group of galaxies with similar mass and size of NGC 1380 gives log₁₀(M_*, ExSitu)=1.063 × log₁₀(M_{*, halo(r <  2 Re)})−0.37 with ℛ = 0.95 and σ = 0.14. With M_{*, halo(r <  2 Re)} = 3.6 ± 0.5 × 10¹⁰ M_⊙, for NGC 1380, we obtained M_*, ExSitu = 10^{10.85 ± 0.15} M_⊙, ($7.1_{-2.0}^{+3.0}\times {10}^{10}$ M_⊙) with a 3σ significant lower limit of 2.5 × 10¹⁰ M_⊙.

Fig. 6. Stellar mass accreted from all galaxies ever merged with NGC 1380 and NGC 1427, i.e., their total ex situ mass. Panel a: same as panel a of Fig. 5, but with stellar mass M*, ExSitu accreted from all satellites on the y axis. The black dots are the subsample of TNG50 galaxies with similar mass and size as NGC 1380. The correlation is stronger, with larger values of the Pearson correlation coefficients, ℛ, and smaller 1σ scatter than the correlation in Fig. 5. The vertical red line with shadows marks the mass and corresponding ±1σ uncertainty of NGC 1380’s hot inner stellar halo. The horizontal dashed line marks the stellar mass accreted from all satellites inferred from the correlations. The thick and thin curves in gray, blue, and green in the background are the running median and ±1σ scatter of the whole sample of TNG50, TNG100, and EAGLE galaxies. Panel b: same as panel a but for NGC 1427.

In a similar fashion, the set of galaxies with similar mass and size of NGC 1427 gives log₁₀(M_*, ExSitu)=1.08 × log₁₀(M_{*, halo(r <  2 Re)})−0.53 with ℛ = 0.83 and σ = 0.22. With M_{*, halo(r <  2 Re)} = (1.4 ± 0.2) × 10¹⁰ M_⊙, for NGC 1427, we inferred M_*, ExSitu = 10^{10.41 ± 0.23} M_⊙, ($2.6_{-1.0}^{1.7}\times {10}^{10}$ M_⊙) with a 3σ significant lower limit of 0.5 × 10¹⁰ M_⊙.

Our results depend insignificantly, within the reported error bars, on the particular subset of simulated galaxies adopted for the inference and in fact do not depend on the specific cosmological galaxy simulation that the relationships have been extracted from either. The correlations fitted for the two subsamples of galaxies chosen for NGC 1380 and NGC 1427 are consistent with each other. In Figs. 5 and 6, we overlap the correlations (Zhu et al. 2022) inferred from all galaxies with 10^10.3 <  M_* <  10^11.6 M_⊙, for three cosmological galaxy simulations whose data are publicly available and that we analyzed in Zhu et al. (2022): TNG50 (Pillepich et al. 2019; Nelson et al. 2019a), TNG100 (Nelson et al. 2019b), and EAGLE (McAlpine et al. 2016). The median as well as the 1σ scatter of these correlations are almost identical in the M_{*, halo(r <  2 Re)} regimes of NGC 1380 and NGC 1427: the difference of M_*, Ex1 and M_*, ExSitu inferred from these three sets of simulations are smaller than 0.1 dex.

4.2. Timing of the past major mergers

The time of the mergers in observed galaxies can be inferred from the stellar age distributions in the different dynamical components. A dynamically cold disk made of stars with λ_z >  0.8 is usually, at least partially, destroyed during a massive merger, resulting in stars of similar age as the surviving disk but on dynamically warmer orbits. Only after the merger, and in the presence of sufficient gas, new stars can form and regrow a stellar disk: this structure can persist for a long time if not disrupted by subsequent massive merger events, whereas minor mergers are expected to only heat up the disk without markedly affecting its circularity distribution (Read et al. 2008). This scenario is supported, for example, by the analysis of data in our Milky Way (Belokurov et al. 2020).

NGC 1380 has a massive dynamically cold disk component. We study the distribution of stars in the circularity λ_z versus stellar age t plane, p(λ_z, t), illustrated by panel a of Fig. 7 for our model of NGC 1380. At any given stellar age t, we can calculate the fraction of disk-like orbits f_disk. Adopting f_disk = 0.5 as the transition from spheroid-dominated to disk-dominated regions, the transition time t|_{fdisk = 0.5} can be associated with the end of the last massive merger. This transition time t|_{fdisk = 0.5} can be well recovered by our model for the galaxies with a massive cold disk at the time of inspection, as is the case for NGC 1380 (see Appendix A). We found such a transition in our model of NGC 1380, as shown in Fig. 7: for ages younger than t ≃ 10 Gyr, stars are suddenly only on cold orbits with λ_z >  0.8 (i.e., f_disk approaches unity, with the transition time of t|_{fdisk = 0.5} = 10 Gyr). We thus inferred that the massive merger responsible for the buildup of NGC 1380’s hot inner stellar halo ended ∼10 Gyr ago. We note that the preexisting cold disk of NGC 1380 was not fully destroyed by the merger, and the surviving disk mainly locates in the inner regions of the galaxy (see also panel b-1 of Fig. 3).

Fig. 7. Timing of the last massive mergers for NGC 1380 and NGC 1427. Panel a: probability distribution of stellar orbits p(λz, t) in the plane of circularity, λz, versus stellar age, t, from our best fitting models of NGC 1380. The fraction of disk orbits, fdisk, is plotted in the top subpanel as a function of t and shows a sharp transition at t ≃ 10 Gyr, indicated by the vertical dashed line. We thus infer that the massive merger responsible for the buildup of NGC 1380’s hot inner stellar halo ended about ∼10 Gyr ago. Panel b: stellar age distribution of the disk, warm, bulge, and hot inner stellar halo components from our best fitting models of NGC 1427. Since the hot inner stellar halo has younger stellar populations than most of the bulge, with the youngest stellar population aged ∼8 Gyr, we infer that the last massive merger responsible for the buildup of the hot inner stellar halo could not have ended before that.

NGC 1427 does not have a massive cold disk today, and so possibly it has not regrown a massive disk after its most massive merger event. Instead, we found that the stars in its hot inner stellar halo are younger than most stars in the bulge; moreover, the hot inner stellar halo holds the youngest stars of the galaxy. According to our understanding of galaxy evolution (see above), inner-halo stars should have formed before the massive merger ended, either in a preexisting disk of the main progenitor or in the merging companion before and during the merger, because of the triggering of star formation. The youngest stars in the hot inner stellar halo therefore set an upper-limit to the end time of the massive merger. Such a scenario is illustrated by a NGC 1427-like galaxy from TNG50 in the appendix. From the stellar age distribution of the hot inner stellar halo of NGC 1427 we obtained (see panel b of Fig. 7) that the massive merger of NGC 1427 has likely ended ≲8 Gyr ago.

NGC 1380 is located in the so-called north-south clump of the Fornax cluster (Iodice et al. 2019b; Nasonova et al. 2011), which comprises a group of galaxies. NGC 1427 is the brightest elliptical galaxy on the eastern side of the Fornax cluster. Both galaxies are thought to be accreted into the potential well of the massive cluster core more than 8 Gyr ago (Iodice et al. 2019b). The galaxy-galaxy merger frequency in cluster environments is expected to be low because of the high relative velocity of their members (Makino & Hut 1997; De Lucia et al. 2004): it is therefore expected that most Fornax galaxies underwent their massive mergers before they fell into the Fornax cluster. The merger of NGC 1380 with its most massive satellite occurring ∼10 Gyr ago is consistent with this scenario. Two other ETGs in the same clump, NGC 1380A (FCC 177) and NGC 1381 (FCC 170), have also been suggested to have accreted a substantial fraction of stars ∼9 Gyr ago (Pinna et al. 2019a,b; Poci et al. 2021). On the other hand, for NGC 1427 and considering that the merger process could last for 1−2 Gyr, the constraint from the stellar age distribution of its hot inner stellar halo does allow the inferred merger event to have started before infall. Whereas we cannot set a strong lower limit on the time of the ending of the most massive merger of NGC 1427, we notice that mergers in a cluster may still happen if two galaxies fell into the cluster within the same group (Tran et al. 2005; Delahaye et al. 2017; Watson et al. 2019).

5. Effects of cluster environment

Since both NGC 1380 and NGC 1427 are in the Fornax cluster, their internal structures might be altered by physical processes like “harassment” or “tidal shocking” in the cluster gravitational potential (e.g., Joshi et al. 2020). TNG50 covers a large variety of group and cluster environments, including ten clusters within the mass range of M₂₀₀ ∼ 10^13.3 − 10^14.3 M_⊙, similar to the Fornax cluster. We can hence directly check with TNG50 whether, and by how much, high-density environments affect the inner stellar halos of satellite galaxies.

In particular, we inspected all the satellite galaxies with M_* >  10⁹ M_⊙, in the second most massive halo of TNG50, which has M₂₀₀ = 9.4 × 10¹³ M_⊙. By comparing the stellar orbit distribution of each galaxy at the time it fell into the cluster to its distribution at z = 0, we find that for ∼80% of the galaxies, the hot inner stellar halo mass M_{*, halo(r <  2 Re)} remains almost unchanged, with differences smaller than 0.1 dex. For the remaining ∼20% of the galaxies, with relatively low mass and that fell into the cluster early, M_{*, halo(r <  2 Re)} decreases significantly as these galaxies get strongly stripped by the cluster. There are a few galaxies whose M_{*, halo(r <  2 Re)} slightly increases, mostly within 0.1 dex. For these systems, their preexisting cold disks were heated up, possibly by harassment, with most of the stars becoming dynamically warm and a small fraction of them becoming hot. There are no satellite galaxies with a significantly hot inner stellar halo component induced after falling into the cluster, due to the fact that no massive mergers happen for these galaxies after falling into the cluster. The massive hot inner stellar halos of the two galaxies studied in this work are thus unlikely to be induced by either harassment or tidal shocking in the cluster environment.

To further illustrate the possible effects of cluster environments on our results, we show the correlations of the hot inner stellar halo mass M_{*, halo(r <  2 Re)} versus accreted stellar masses M_*, Ex1 and M_*, ExSitu for the central and satellite galaxies from TNG50 separately. We use all TNG50 galaxies with M_* = 10^10.3 − 10^11.6 M_⊙, except 37 with ongoing mergers – the same sample as used in Zhu et al. (2022). We quantify the stripped mass of a galaxy by comparing its present-day stellar mass M_{*, z = 0} with the maximum stellar mass it has ever reached M_∗,max. All satellite galaxies with M_∗,z=0 < 0.8 M_∗,max are considered to be strongly stripped. In Fig. 8 we can see that the correlations are almost the same for central and satellite galaxies except the strongly stripped ones. Most of the strongly stripped galaxies have M_{*, halo(r <  2 Re)} <  10¹⁰ M_⊙, they tend to be located below the median of the correlations, and some of them are located significantly below the lower “−1σ” curves.

Fig. 8. Correlations of hot inner stellar halo mass, M*, halo(r <  2 Re), versus accreted stellar mass M*, Ex1 (left) and M*, ExSitu (right) for central and satellite galaxies separately. All TNG50 galaxies with M* = 1010.3 − 1011.6 M⊙, are included. The black dots are central galaxies, and the red dots are satellites. The thick black and red curves are the running median for central and satellite galaxies, respectively, while the thin dashed curves are the corresponding ±1σ scatter. The satellites with more than 20% of stellar mass stripped are marked by magenta crosses. The correlations are almost the same for central and satellite galaxies. The strongly stripped satellites tend to be located below the median of the correlations, and some of them deviate from the correlations significantly below the lower “−1σ” curves.

Instead of M_*, Ex1 and M_*, ExSitu, we can quantify a merger by the secondary progenitor mass defined as the maximum of the stellar mass it has ever reached before the merger ended (Rodriguez-Gomez et al. 2016). We use M_*, Sat1 to indicate the stellar mass of the most massive secondary progenitor a galaxy has ever accreted, and M_*, allSat to indicate the sum of all the secondary progenitors more massive than 1% stellar mass of the current galaxy. M_*, Ex1 is highly correlated and larger than M_*, Sat1 by a factor of ∼30% for all the central galaxies and for all the satellite galaxies that are not strongly stripped³. On the other hand, for the strongly stripped galaxies, M_*, Ex1 calculated from the current galaxy are significantly smaller than the original accreted stellar mass, and also smaller than M_*, Sat1. In Fig. 9 we show the correlations by using M_∗,allSat and M_∗,Sat1 instead of M_*, ExSitu and M_*, Ex1. The strongly stripped galaxies, which were located below the −1σ curves in Fig. 8, move upward, and tend to follow the relations in Fig. 9. Some of the strongly stripped galaxies become outliers above the upper “+1σ” curves.

Fig. 9. Correlations from using the secondary progenitor masses M∗,allSat and M∗,Sat1 instead of the ex situ stellar masses in present-day galaxies, M∗,ExSitu and M∗,Ex1. Symbols are the same as in Fig. 8. The strongly stripped galaxies move upward, and some of them become outliers above the upper “+1σ” curves.

In summary, for all central and most satellite galaxies, M_*, Ex1 and M_*, Sat1 are equally valid to quantify their past merger events, while M_{*, halo(r <  2 Re)} defined in the current galaxy is a good indicator for the merger mass. However, if a satellite galaxy is strongly stripped, then both M_*, Ex1 and M_{*, halo(r <  2 Re)} can be affected, and the former more than the latter. The accreted stellar mass M_*, Ex1 inferred from the median correlation in Fig. 8 could be higher than M_*, Ex1 it has in the current galaxy, but lower than the original M_*, Ex1 it once had. Similarly, M_*, Sat1 inferred from the median correlation in Fig. 9 is likely to be underestimated if the hot inner stellar halo was strongly stripped.

For our two galaxies in the Fornax cluster, NGC 1380 is a massive lenticular galaxy with M_* ∼ 1.8 × 10¹¹ M_⊙, with median size for its stellar mass, and with an extended disk. Its outer isophotoes at R >  3 arcmin appear twisted with respect to the inner, diskier regions (Iodice et al. 2019a), which could be the result of the major merger. All these features suggest that NGC 1380 is unlikely to be strongly stripped and the results presented above should not be affected by environmental effects. NGC 1427 is a massive elliptical galaxy with M_* ∼ 5 × 10¹⁰ M_⊙, also with median size for its stellar mass, and with a massive hot inner stellar halos with M_{*, halo(r <  2 Re)} >  10¹⁰ M_⊙. The outer isophotes of NGC 1427 appear asymmetric, suggesting it could be partly stripped. However, its inner regions – and thus also the hot inner stellar halo – are unlikely to be significantly affected. We expect that M_{*, halo(r <  2 Re)} is still a good indicator for its past merger mass, although the ex situ stellar mass in the current galaxy might actually be lower.

6. Summary

We have discovered that the two ETGs in the Fornax cluster, NGC 1380 and NGC 1427, both underwent an ancient massive merger. We have constrained the merger mass and time.

By applying our recently developed population-orbit superposition method to the observed luminosity density, stellar kinematic, age, and metallicity maps, we obtained the stellar orbits, the age, and the metallicity distribution of each galaxy. We then decomposed each galaxy into four components, which included a dynamically hot inner stellar halo as a relic of past massive mergers. By comparing to analogs from cosmological galaxy simulations, chiefly TNG50, we used the robustly measured mass of the hot inner stellar halo to infer the stellar mass from the ever-accreted most massive satellite, M_*, Ex1, as well as the total ex situ stellar mass, M_*, ExSitu, for these two galaxies. We further found that the last merger of NGC 1380 ended ∼10 Gyr ago based on the stellar age distribution of its regrown dynamically cold disk, in a way comparable to how the last massive merger time of the Milky Way could be inferred (Belokurov et al. 2020). The merger in NGC 1427 ended t ≲ 8 Gyr ago based on the stellar populations in its hot inner stellar halo, in a way comparable to how the merger time in M 31 was estimated (D’Souza & Bell 2018b). While the quantification of the satellite masses and merger times in both the Milky Way and M 31 rely on observations of individual stars, we overcame these limitations by applying our new population-orbit superposition method to IFU data.

In Table 1 we list the mass and time of the most massive merger events of NGC 1380 and NGC 1427 in comparison with those of the Milky Way and M 31. The Milky Way, M 31, and NGC 1380 are the only three galaxies for which both the merger mass and time have been quantitatively uncovered so far. The mass of the most massive merger of NGC 1427 is comparable to that of M 31, while the mass of the most massive merger of NGC 1380 is twice as large. However, the mergers in NGC 1380 and NGC 1427 happened much earlier than in M 31. Given the subsequent evolution in the cluster environment, the high encounter velocities and long relaxation times have resulted in smoother structural distributions in NGC 1380 and NGC 1427 than in M 31. The most massive merger in the Milky Way is similarly old but is smaller, with a stellar mass of 3−6 × 10⁸ M_⊙, which is only about ∼1.5 per cent of NGC 1380’s and ∼3 per cent of NGC 1427’s progenitor satellite stellar mass. The merger event in NGC 1380 is the oldest and most massive one uncovered in nearby galaxies so far.

Integral-field-unit data similar to what we have used for NGC 1380 and NGC 1427 could be obtained for hundreds of nearby galaxies through observations with instruments such as MUSE. Our discovery thus opens a new way to quantify the merger history of nearby galaxies and places important constraints on the assembly of galaxies as a function of environment and in a cosmological context.

Acknowledgments

We are grateful for discussions with Illustris TNG team members Dylan Nelson, Lars Hernquist, Ruediger Pakmor and Mark Vogelsberger. L.Z. acknowledges the support from the National Key R&D Program of China under grant No. 2018YFA0404501 and National Natural Science Foundation of China under grant No. Y945271001. G.v.d.V. 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). E.M.C. is supported by MIUR grant PRIN 2017 20173ML3WW_ and Padua University grants DOR1885254/18, DOR1935272/19, and DOR2013080/20. F.P., I.M.N., and J.F.B. 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. TNG50 was realised with compute time granted by the Gauss Centre for Super-computing (GCS), under the GCS Large-Scale Project GCS-DWAR (2016; PIs Nelson/Pillepich).

References

Appendix A: Method validation

We previously tested the ability of our population-orbit method to recover the age and metallicity properties of different orbital components (Zhu et al. 2020). However, there we did so only for simulated spiral galaxies that do not have without a significant hot inner stellar halo within the data coverage. In this Appendix, we extend such tests.

The major degeneracy that our model could present is between the disk and the halo, as a face-on disk might look similar to a halo in the kinematic maps. NGC 1427 has a small disk in the inner regions with r <  10″. The disk fraction in the outer regions is very small. Therefore, the disk-halo degeneracy should not be a problem for NGC 1427. For NGC 1380, on the other hand, we performed a series of tests using mock data from four simulated galaxies from TNG50: TNG50 468590, TNG50 375073, TNG50 2, and TNG50 220596, which we chose due to the similarity between their hot inner stellar halo mass and disk mass and those of NGC 1380.

The mock data were created in a similar way as done before (Zhu et al. 2020), but restricted to have the same inclination angle, ϑ = 77°, and similar data coverage and data quality as for NGC 1380. We first projected a simulated galaxies into the observational plane with inclination angle ϑ = 77°, placed it at a distance of d = 25 Mpc, and observed it with a pixel size of 1″. Then we summed up the mass of the stellar particles within each pixel to obtain a surface mass density map (i.e., a mock image of the galaxy) by assuming that light equals mass.

We chose the area with −40″ <  x <  40″ and −35″ <  y <  35″ to create the kinematic maps. This data coverage is chosen to be similar to the observed kinematic maps of NGC 1380, where x is along the major axis and y is along the minor axis of the galaxy. According to the number of stellar particles in each pixel, we perform a Voronoi binning process to reach a target (S/N)_T = 50, assuming Poission noise $\sim \sqrt{N_{\mathrm{particles}}}$. Given all the particles in each Voronoi bin, we derive the mass-weighted mean velocity, velocity dispersion, h₃, and h₄ by fitting a GH function (van der Marel & Franx 1993) to the stellar line-of-sight velocity distribution (LOSVD), as well as we calculate mass-weighted average age (t) and metallicity (Z/Z_sun). We note that here the mean velocity and velocity dispersion are not directly calculated from the LOSVD, but they are actually the parameters of the base Gaussian in the GH function.

After the Voronoi binning, the spatial resolution of our mock data is ∼150−1000 pc, which is comparable to that of the kinematic data (binned with target (S/N)_T = 100) from the Fornax 3D project (Sarzi et al. 2018). We use a simple function inferred from the CALIFA data to construct the errors for kinematic maps (Tsatsi et al. 2015), with the errors proportional to $\frac{{(S/N)}_T}{{(S/N)}_{\mathrm{bin}}}(1+1.4log{N}_{\mathrm{pix}})$, where ${(S/N)}_{\mathrm{bin}}=\sqrt{N_{\mathrm{particles}}}$, and N_pix is the number of pixels in a Voronoi bin. We have larger errors for kinematics in the area with lower surface density, where more pixels are included in each bin in order to reach the target (S/N)_T. For age and metallicity, the tests on the full-spectrum fitting to mock spectra of (S/N)_T = 40 return random errors of 10% (Pinna et al. 2019a). The errors could be lower for spectra with higher S/N, while it could be higher for real spectra due to possible systematic effects. As proof-of-concept, we adopt relative errors of 10% for age and metallicity. So, to include the noise, the kinematics, age and metallicity maps are perturbed by random numbers, normally distributed with dispersions equal to the observational errors. We find that the error maps of the mock data are similar to the data of MUSE observations for galaxies in the Fornax 3D project (Sarzi et al. 2018).

We take each mock data set as a real observed galaxy and apply the population-orbit superposition model to it. In Fig. A.1 we show the mock data created from TNG50 468590 and our best fitting population-orbit model to it. The model fits well the data in great detail.

Fig. A.1. Population-orbit model of the simulated galaxy TNG50 468590. The first column shows the mock data with, from top to bottom, surface mass density, mean velocity, velocity dispersion, GH coefficients h3, h4, mass-weighted age, and metallicity maps. The second column is the best fitting population-orbit model, and the third column shows the residuals of the data minus model, divided by the errors assigned assuming similar quality data as for NGC 1380.

A.1. Recovery of the four components

In Fig. A.2 we show that the probability density distribution p(λ_z, r), age distribution t(λ_z, r), and metallicity distribution Z(λ_z, r) directly measured from the simulations are well matched by those extracted from the best fitting population-orbit models. There is a discrepancy in the inner parts, where the bulge in the simulated galaxies is dominated by box orbits with λ_z = 0, whereas part of the box orbits in the model are actually represented by tube orbits with nonzero but opposite values of λ_z. This discrepancy, however, does not affect our division into four components: disk (λ_z >  0.8), bulge (λ_z <  0.8, r <  r_cut), warm component (0.5 <  λ_z <  0.8, r >  r_cut), and hot inner stellar halo (λ_z <  0.5, r_cut <  r <  2R_e), where r_cut = 3.5 kpc is chosen for all cases. The mass fraction, average age and metallicity of the disk, bulge, warm component, and hot inner stellar halo within a 2R_e radius for the four simulated galaxies are given in Table A.1. Our model excellently recovers the mass fraction, age and metallicity of the hot inner stellar halo, even though the statistical errors obtained from the 1σ scatter of the best fitting models are too small to cover the true values in some cases. The absolute difference between the true values and our model recovered mean values are still small.

Fig. A.2. Verification of the population-orbit method with mock data created from galaxies of the TNG50 simulation: TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2 from Panel (a) to Panel (d). For each galaxy, the left three plots show the distributions directly from the simulation, which are well matched by those extracted from the best fitting population-orbit models in the right three plots. From left to right, we show the distributions in phase space of the time-averaged radius, r, versus circularity, λz, of orbital probability density p(r, λz), age t(r, λz), and metallicity Z(λz, r), all shown only within 2 Re. The dashed lines indicate our orbital-based division into four components: disk, bulge, warm component, and hot inner stellar halo. The corresponding mass fractions within a 2 Re radius are listed in Table A.1.

Fig. A.3. Verification of the bulge versus hot inner stellar halo decomposition. Shown are the probability density phot (top panel) and metallicity Zhot/Z⊙ (bottom panel) of the hot orbits as a function of radius for TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2, from Panel (a) to Panel (d). In each panel, the red curve is the true profile from the simulation, and the thick and dashed black curves are the mean and 1σ scatter of our best fitting models. The vertical lines mark the position of r = 3.5 kpc, where we separate the bulge from hot inner stellar halo.

Fig. A.4. Verification of the method to obtain merger time using the circularity-age distribution. Shown are probability density distributions in phase space of stellar age, t, versus circularity, λz, for TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2, from Panel (a) to Panel (d). In each panel, the left plot is the true distribution for the stellar orbits within 2 Re directly obtained from the simulation, and the vertical dashed line marks the merger end time of the last massive merger the galaxy has experienced. The merger end time is consistent with the time t|fdisk = 0.5 in TNG50 468590, TNG50 375073, and TNG50 220596, with still significant star formation after the last massive merger, and it is consistent with the end of star formation in TNG50 2. The right plot is the distribution from our best fitting model, and the vertical dashed line marks the time t|fdisk = 0.5 in the model for TNG50 468590, TNG50 375073, and TNG50 220596, while it marks the end of star formation in TNG50 2.

In Fig. A.3 we display the probability density distribution p_hot and average metallicity Z_hot as a function of radius r by summing all the hot orbits with λ_z <  0.5. The red curves are directly calculated from the simulations, while the thick and dashed black curves are the mean and 1σ scatter from our best fitting models. The model generally matches the radial profiles for both the probability density and metallicity. Both p_hot and Z_hot sharply decrease from an inner peak, and flatten at larger radii.

From the above analysis, we find that, although our data have limited extension along the galaxy minor axis, we still are able to uncover the intrinsic stellar orbit distributions quite well. This is thanks to the high data quality that provide reliable high-order GH coefficients h₃ and h₄, which in turn encode the information of the combination of disk and hot inner stellar halo. The age and metallicity of stellar orbits in the hot inner stellar halo from our model still have relatively large fluctuations. Having data extending more along the galaxy minor axis should improve the constraints on the stellar age and metallicity distribution of the hot inner stellar halo.

Fig. A.5. Verification of the method to assign a time limit to the merger from the age distributions of the galaxy components. Shown are the stellar age distributions of the disk, warm, bulge, and hot inner stellar halo components for TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2 from Panel (a) to Panel (d). In each panel, the left plot is the true stellar age distribution from the stellar particles in the simulation, the right plot is the distribution from our best fitting model, and the vertical dashed line marks the end of the last massive merger.

A.2. Recovery of the merger’s end time

We further check how well we can recover the correlation between stellar age t and circularity λ_z. In Fig. A.4 we show the probability density of orbits in the circularity λ_z versus stellar age t, p(λ_z, t), of the four simulated galaxies elected for the testing (panels a–d). In each panel, the left plot gives the true distribution from the simulation and the right plot shows our best fitting model. We know every aspect of the merger history of the simulated galaxies: the real end time of the last massive merger of each simulated galaxy is marked in Fig. A.4. In three of the four galaxies, TNG50 468590, TNG50 375073, and TNG50 220596, there is substantial star formation after the last massive merger. Stars formed after the last massive mergers rarely ended on non-disk orbits. In all these cases, a sharp transition of stars from non-disk dominated orbits to disk dominated orbits can be clearly identified: therefore we can confidently associate this transition with the last massive merger end time. The time of t|_{fdisk = 0.5} is a good proxy for the merger end time.

In practice, t|_{fdisk = 0.5} is well recovered by our model in TNG50 468590 and TNG50 375073, which have massive cold disks formed after the last massive merger, although it is not as sharp as in the simulation data due to uncertainties of the stellar age of orbits in our model. In TNG50 220596, only a small fraction of disk formed after the merger: even though our model generally matches the distribution p(λ_z, t), the transition time t|_{fdisk = 0.5} is hard to recover. In TNG50 2, there is no disk formed after the last merger: throughout the stellar age range, the cold disk never becomes dominating. The merger end time is consistent with the star formation end time. Our model generally recovers the distribution p(λ_z, t), and the star formation end time. In both TNG50 220596 and TNG50 2, a massive merger occurring at low redshift partially heated up the previous cold disk and produced a significant warm component. The stars in such a warm component exhibit a wide range of stellar ages, similar to that of the remnant cold orbits (Fig. A.4). However, they are on average younger than the stars on dynamical hot orbits with λ_z ∼ 0. While NGC 1380 has a massive cold disk, which is younger than any other component, the warm component in NGC 1380 is as old as the stars on hot orbits, like those in TNG50 468590 and 375073.

In Fig. A.5 we analyze the recovery of age distribution of the four components. In the left plot of each galaxy, we show the stellar age distribution of particles belonging to each of the four components in the simulation. In the right plot, we show the stellar age distribution of orbital bundles belonging to each component from our best fitting models. Our model is able to roughly recover the true age distribution of each component.

In these four simulated galaxies, inner-halo stars are old and disk stars are young: in some cases, for example TNG50 468590 and TNG50 375093, our model introduces some young stars in the inner-halo due to the degeneracy with disk. Similarly for NGC 1380, we get a small fraction of inner-halo stars with t <  10 Gyr, which are likely induced by such degeneracy. NGC 1380 is thus still consistent with the scenario that the inner-halo stars are formed before the merger end time of 10 Gyr.

In these four simulated galaxies, disk is the youngest component, the youngest stars in the disk being 0, 0, 2, and 2 Gyr old in TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2, respectively: from the best fitting model, we obtain 2, 0, 4, and 3 Gyr. The declination of probability density toward younger ages is sharper in our model than in the real distribution from the simulation. Similarly for the hot inner stellar halo in NGC 1427, which is the youngest component in that galaxy. The star formation end time of 8 Gyr obtained by our model is likely an upper limit, and the real star formation end time could be 1–2 Gyr younger. So the merger end time of t ≲ 8 Gyr is a relaxed upper limit.

Appendix B: Uncertainties on the stellar mass of the hot inner stellar halo

We have shown that the mass fraction of the hot inner stellar halo is well recovered by our population-orbit superposition model fitted to the mock data of NGC 1380-like simulated galaxies from TNG50. The mass fractions of the hot inner stellar halo f_halo of the four TNG50 test galaxies is 0.333, 0.221, 0.214, and 0.30, respectively: we get 0.351 ± 0.004, 0.228 ± 0.005, 0.244 ± 0.008 and 0.31 ± 0.01 from the best fitting models. Here the error is the statistical error taken 1σ scatter of the best fitting models. The statistical error is very small, and the true value is not always covered within the 1σ statistical error. As the absolute difference between the true and model-estimated values is 0.018, 0.007, 0.01, and 0.03 for the four simulated galaxies, we consider the average of these four values as a systematic error, which gives 0.016. For NGC 1380, we obtained the luminosity fraction of hot inner stellar halo f_halo within 2R_e as 0.267 ± 0.006: considering a systematic error of 0.016, we have f_halo = 0.267 ± 0.017. The error of 0.017 is determined by $d=\sqrt{d_1^2+d_2^2}$, where d₁ = 0.006 is the statistical error and d₂ = 0.016 is the systematic error. For NGC 1427, we obtained the luminosity fraction of hot inner stellar halo within 2R_e as 0.34 ± 0.02, and considering a systematic error of 0.016, we have f_halo = 0.34 ± 0.03.

Whereas orbit probabilities in the simulated galaxies are mass-weighted, those for NGC 1380 and NGC 1427 are r-band luminosity weighted. For NGC 1380, under the simplest assumption of best fitting constant stellar mass-to-light ratio M_*/L_r = 2.7 ± 0.3 M_⊙/L_⊙, from our dynamical models, we estimate M_{*, halo(r <  2 Re)} = L_r(r <  2 R_e)×M_*/L_r × f_halo = (3.6 ± 0.5) × 10¹⁰ M_⊙. The uncertainty is calculated by the error propagation formula considering the uncertainties in M_*/L_r and f_halo. Similarly for NGC 1427, we arrive at M_{*, halo(r <  2 Re)} = (1.4 ± 0.2) × 10¹⁰ M_⊙, with M_*/L_r = 2.5 ± 0.2 given by our dynamical model. M_{*, halo(r <  2 Re)} will be used throughout to infer the masses of the mergers.

In NGC 1380, we know from a previous study on stellar population synthesis that M_*/L_r varies from 4.8 ± 0.7 M_⊙/L_⊙, in the bulge-dominated inner regions to 3.5 ± 0.5 M_⊙/L_⊙, in the inner-stellar-halo-dominated outer regions by considering a variable initial mass function (Martín-Navarro et al. 2019). The adopted error is the 1σ scatter of M/L_r in the corresponding regions. We tested in previous work (Zhu et al. 2018a) that assuming a constant mass-to-light ratio (even if it actually varies) to infer the gravitational potential does not affect the stellar orbit distribution in the dynamical model (Zhu et al. 2018a; Poci et al. 2019), and thus does not affect the luminosity fraction of hot inner stellar halo f_halo. Furthermore, the M_*/L_r ratios from stellar population synthesis in the halo-dominated outer regions are roughly consistent with the M_*/L_r ratios inferred from our dynamical model. With M_*/L_r = 3.5 ± 0.5 M_⊙/L_⊙, we have M_{*, halo(r <  2 Re)} = (4.4 ± 0.8) × 10¹⁰ M_⊙, which is still within 1σ of our quoted reference value.

Appendix C: A NGC 1427-like galaxy in TNG50

Fig. C.1. Probability density distributions in phase space of stellar age, t, versus circularity, λz, for stars with different origins in TNG50 465255. From left to right, they are in situ stars formed before the merger ended, stars accreted from the secondary progenitor, and stars formed in situ after the merger ended. The dashed black and blue lines indicate the merger starting and ending time.

Fig. C.2. Some internal properties of TNG50 465255. Panel (a): Probability density of stars on dynamical hot orbits with λz <  0.5 as a function of radius, Panel (b): Metallicity distribution of stars on the hot orbits as a function of radius. Panel (c): Stellar age distribution of different components. Stars in the hot inner stellar halo are younger than those in the bulge, similar to what we found for NGC 1427.

NGC 1427 does not have a significant cold disk and its hot inner stellar halo is younger than the bulge, which is different from NGC 1380. While the galaxies we chose for case studies in Zhu et al. (2022) are all NGC 1380 analogs, here we present a NGC 1427-like galaxy taken from TNG50. About 25% of TNG50 galaxies have hot inner stellar halos younger than bulges. One of these is TNG50 465255 with M_* = 9.3 × 10¹⁰ M_⊙, R_e = 4 kpc, and a hot inner stellar halo mass of M_{*, halo(r <  2 Re)} = 1.2 × 10¹⁰ M_⊙. The galaxy has experienced a massive merger started at t ∼ 6 Gyr and ended at t ∼ 4.8 Gyr with a stellar mass accreted from the secondary progenitor of M_∗,Ex1 = 1.7 × 10¹⁰ M_⊙. No other merger more massive than 0.1 M_∗,Ex1 happened after that. The total ex situ stellar mass is M_∗,ExSitu = 2.2 × 10¹⁰ M_⊙.

We split the stars in the galaxy at z = 0 into three categories according to their origins: (1) stars formed in situ before the merger ended, (2) stars accreted from the secondary progenitor, and (3) stars formed in situ after the merger ended. In Fig. C.1 we show the probability density distributions in phase space of stellar age t versus circularity λ_z for each part. The in situ stars formed before the merger ended have a wide range of λ_z in the galaxy at z = 0. The dynamical hot stars contribute mostly to the bulge and they are very old. The stars accreted from the secondary progenitor are younger and mostly contribute to the hot inner stellar halo. After the merger ended, only very few stars formed on dynamically cold orbits with little gas left after the merger. At the end, TNG50 465255 is an elliptical galaxy.

We separate the galaxy into four components: cold disk, warm component, bulge and hot inner stellar halo in the same way as we did for NGC 1427. In the left panels of Fig. C.2, we show the probability density and metallicity distribution of the hot orbits (with λ_z <  0.5) as a function of radius. Here we use radial bin size of ∼1 kpc to avoid the dip of probability density at the inner-most bin. We separate the bulge and hot inner stellar halo at r = 3.5 kpc as we did for all the galaxies (Zhu et al. 2022), although the transition in density is not sharp and the metallicity has already flatten at smaller radius, similarly to NGC 1427. The stellar age distribution of the four components is shown in the right panel of Fig. C.2. Similar to NGC 1427, TNG50 465255 has a younger inner hot stellar halo than the bulge. For this galaxy, the age of the youngest star in the hot inner stellar halo indicates the merger ending time.

All Tables

All Figures

Fig. 1. Population-orbit model of NGC 1380. Panel a: deep r-band image of NGC 1380 from the FDS with the VLT Survey Telescope (Iodice et al. 2019a). The three squares indicate the three pointings, each with a 1′ × 1′ square field-of-view, of spectroscopic observations with the VLT MUSE instrument. The black dashed ellipse corresponds to the isophote at μB = 25 mag arcsec−2. Panel b: data from the VLT/MUSE observations (first column), the best fitting population-orbit model (second column), and the residuals of the data minus model divided by the observational errors (third column). From top to bottom: surface brightness, mean velocity, velocity dispersion, GH coefficients h3 and h4, and the light-weighted age and metallicity maps of the stars in NGC 1380. We note that kinematic maps were derived from spectra binned with a different scheme than the age and metallicity maps.

In the text

	Fig. 2. Population-orbit model of NGC 1427. See the caption of Fig. 1 for a description.
In the text

Fig. 3. Orbital decompositions of NGC 1380 and NGC 1427. Panel a-1: probability density distribution of the stellar orbits p(r, λz) in the phase space of time-averaged radius r versus circularity λz in units of stellar mass per unit area in the phase space. We show the model within 140 arcsecs (or ∼14 kpc) for NGC 1380. The stellar masses are normalized such that they sum to unity within the data coverage. Panels b-1 and c-1: stellar age t(r, λz) and metallicity Z(t, λz) distribution of the orbits in the same phase space. All distributions p(r, λz), t(r, λz), and Z(t, λz) are averages of multiple best fitting models that fall within the 1σ confidence level of the model hyperparameters. The dashed lines indicate our orbit-based division into four components: disk, warm, bulge, and hot inner stellar halo. Panels a-2, b-2, c-2: same but for NGC 1427, which does not have an extended disk.

In the text

Fig. 4. Probability density and metallicity profiles. Panel a-1: probability density of dynamically hot orbits with λz <  0.5 (red) and warm orbits with 0.5 <  λz <  0.8 (orange) as a function of radius for NGC 1380. The red and orange thick and dashed curves are the mean and 1σ scatter of our best fitting models. The vertical line at r = 3.5 kpc marks the separation of the bulge from the hot inner stellar halo. Panel b-1: metallicity distribution of hot and warm orbits as a function of radius for NGC 1380. Panels a-2 and b-2: same but for NGC 1427.

In the text

Fig. 5. Stellar masses of the most massive satellites (now-destroyed) that merged with NGC 1380 and NGC 1427. Panel a: correlation from cosmological simulations between the mass of the hot inner stellar halo, M*, halo(r <  2 Re), and the stellar mass, M*, Ex1, of the most massive merged satellite. The black dots are a subsample of TNG50 galaxies with similar mass and size as NGC 1380, the thick black line is a linear fit to the points, and the thin black lines indicate the ±1σ scatter. The Pearson correlation coefficients, ℛ, between M*, halo(< 2 Re) and M*, Ex1 and 1σ scatter against the linear fit for the TNG50 subsample are indicated. The vertical red line with shadows marks the mass and corresponding ±1σ uncertainty of NGC 1380’s hot inner stellar halo. The horizontal dashed line marks the satellite(s) mass inferred from the correlations. The thick and thin curves in gray, blue, and green in the background are the running median and ±1σ scatter of the whole sample of TNG50, TNG100, and EAGLE galaxies. Panel b: same to panel a, but for the subsample of TNG50 galaxies with similar mass and size to NGC 1427. The vertical red line marks the mass of NGC 1427’s hot inner stellar halo, and the horizontal dashed line marks the mass of the satellite(s) inferred from the correlations.

In the text

Fig. 6. Stellar mass accreted from all galaxies ever merged with NGC 1380 and NGC 1427, i.e., their total ex situ mass. Panel a: same as panel a of Fig. 5, but with stellar mass M*, ExSitu accreted from all satellites on the y axis. The black dots are the subsample of TNG50 galaxies with similar mass and size as NGC 1380. The correlation is stronger, with larger values of the Pearson correlation coefficients, ℛ, and smaller 1σ scatter than the correlation in Fig. 5. The vertical red line with shadows marks the mass and corresponding ±1σ uncertainty of NGC 1380’s hot inner stellar halo. The horizontal dashed line marks the stellar mass accreted from all satellites inferred from the correlations. The thick and thin curves in gray, blue, and green in the background are the running median and ±1σ scatter of the whole sample of TNG50, TNG100, and EAGLE galaxies. Panel b: same as panel a but for NGC 1427.

In the text

Fig. 7. Timing of the last massive mergers for NGC 1380 and NGC 1427. Panel a: probability distribution of stellar orbits p(λz, t) in the plane of circularity, λz, versus stellar age, t, from our best fitting models of NGC 1380. The fraction of disk orbits, fdisk, is plotted in the top subpanel as a function of t and shows a sharp transition at t ≃ 10 Gyr, indicated by the vertical dashed line. We thus infer that the massive merger responsible for the buildup of NGC 1380’s hot inner stellar halo ended about ∼10 Gyr ago. Panel b: stellar age distribution of the disk, warm, bulge, and hot inner stellar halo components from our best fitting models of NGC 1427. Since the hot inner stellar halo has younger stellar populations than most of the bulge, with the youngest stellar population aged ∼8 Gyr, we infer that the last massive merger responsible for the buildup of the hot inner stellar halo could not have ended before that.

In the text

Fig. 8. Correlations of hot inner stellar halo mass, M*, halo(r <  2 Re), versus accreted stellar mass M*, Ex1 (left) and M*, ExSitu (right) for central and satellite galaxies separately. All TNG50 galaxies with M* = 1010.3 − 1011.6 M⊙, are included. The black dots are central galaxies, and the red dots are satellites. The thick black and red curves are the running median for central and satellite galaxies, respectively, while the thin dashed curves are the corresponding ±1σ scatter. The satellites with more than 20% of stellar mass stripped are marked by magenta crosses. The correlations are almost the same for central and satellite galaxies. The strongly stripped satellites tend to be located below the median of the correlations, and some of them deviate from the correlations significantly below the lower “−1σ” curves.

In the text

	Fig. 9. Correlations from using the secondary progenitor masses M∗,allSat and M∗,Sat1 instead of the ex situ stellar masses in present-day galaxies, M∗,ExSitu and M∗,Ex1. Symbols are the same as in Fig. 8. The strongly stripped galaxies move upward, and some of them become outliers above the upper “+1σ” curves.
In the text

Fig. A.1. Population-orbit model of the simulated galaxy TNG50 468590. The first column shows the mock data with, from top to bottom, surface mass density, mean velocity, velocity dispersion, GH coefficients h3, h4, mass-weighted age, and metallicity maps. The second column is the best fitting population-orbit model, and the third column shows the residuals of the data minus model, divided by the errors assigned assuming similar quality data as for NGC 1380.

In the text

Fig. A.2. Verification of the population-orbit method with mock data created from galaxies of the TNG50 simulation: TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2 from Panel (a) to Panel (d). For each galaxy, the left three plots show the distributions directly from the simulation, which are well matched by those extracted from the best fitting population-orbit models in the right three plots. From left to right, we show the distributions in phase space of the time-averaged radius, r, versus circularity, λz, of orbital probability density p(r, λz), age t(r, λz), and metallicity Z(λz, r), all shown only within 2 Re. The dashed lines indicate our orbital-based division into four components: disk, bulge, warm component, and hot inner stellar halo. The corresponding mass fractions within a 2 Re radius are listed in Table A.1.

In the text

Fig. A.3. Verification of the bulge versus hot inner stellar halo decomposition. Shown are the probability density phot (top panel) and metallicity Zhot/Z⊙ (bottom panel) of the hot orbits as a function of radius for TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2, from Panel (a) to Panel (d). In each panel, the red curve is the true profile from the simulation, and the thick and dashed black curves are the mean and 1σ scatter of our best fitting models. The vertical lines mark the position of r = 3.5 kpc, where we separate the bulge from hot inner stellar halo.

In the text

Fig. A.4. Verification of the method to obtain merger time using the circularity-age distribution. Shown are probability density distributions in phase space of stellar age, t, versus circularity, λz, for TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2, from Panel (a) to Panel (d). In each panel, the left plot is the true distribution for the stellar orbits within 2 Re directly obtained from the simulation, and the vertical dashed line marks the merger end time of the last massive merger the galaxy has experienced. The merger end time is consistent with the time t|fdisk = 0.5 in TNG50 468590, TNG50 375073, and TNG50 220596, with still significant star formation after the last massive merger, and it is consistent with the end of star formation in TNG50 2. The right plot is the distribution from our best fitting model, and the vertical dashed line marks the time t|fdisk = 0.5 in the model for TNG50 468590, TNG50 375073, and TNG50 220596, while it marks the end of star formation in TNG50 2.

In the text

Fig. A.5. Verification of the method to assign a time limit to the merger from the age distributions of the galaxy components. Shown are the stellar age distributions of the disk, warm, bulge, and hot inner stellar halo components for TNG50 468590, TNG50 375073, TNG50 220596, and TNG50 2 from Panel (a) to Panel (d). In each panel, the left plot is the true stellar age distribution from the stellar particles in the simulation, the right plot is the distribution from our best fitting model, and the vertical dashed line marks the end of the last massive merger.

In the text

	Fig. C.1. Probability density distributions in phase space of stellar age, t, versus circularity, λz, for stars with different origins in TNG50 465255. From left to right, they are in situ stars formed before the merger ended, stars accreted from the secondary progenitor, and stars formed in situ after the merger ended. The dashed black and blue lines indicate the merger starting and ending time.
In the text

Fig. C.2. Some internal properties of TNG50 465255. Panel (a): Probability density of stars on dynamical hot orbits with λz <  0.5 as a function of radius, Panel (b): Metallicity distribution of stars on the hot orbits as a function of radius. Panel (c): Stellar age distribution of different components. Stars in the hot inner stellar halo are younger than those in the bulge, similar to what we found for NGC 1427.

In the text