© ESO 2022

1. Introduction

Galaxies can experience a diversity of assembly and merger histories, which in turn are thought to lead to the observed variety in galaxy star formation rates, stellar morphologies, and kinematics. Galaxies with quiescent merger histories are thought to be likely dominated by fast rotating, flat stellar disks with recent star formation (e.g., Fall & Efstathiou 1980). On the other hand, major (non-gas-rich) mergers can play a key role in destroying stellar disks and in creating dynamically hot spheroidal components in the remnant galaxies (Cox et al. 2006; Hoffman et al. 2010; Bois et al. 2010, 2011; Naab et al. 2014; Pillepich et al. 2015). Minor mergers are thought to be efficient in building the outer stellar halos (Naab et al. 2009).

With deep photometric imaging, observed galaxies can be decomposed into multiple stellar components, such as a disk, bulge, bar, and halo. Great efforts over the past few years have been made to find the connection between galactic substructures and merger events. For example, it has been suggested that the bulges of Milky Way-mass galaxies have not been built up by mergers (Kormendy & Kennicutt 2004), or at least not exclusively (Bell et al. 2017). On the other hand, deep imaging does not yet seem to provide enough information to infer the details of the accreted satellites that formed the stellar halos (Spavone et al. 2020), and their very separation from other galactic components may be affected by large uncertainties in comparison to kinematically based decomposition methods (Du et al. 2020).

Over the last two decades, the rapidly growing field of integral-field-unit (IFU) spectroscopy has significantly deepened our knowledge of stellar structures in galaxies across a wide range of masses and Hubble types. A number of IFU surveys, such as SAURON (Davies et al. 2001), ATLAS-3D (Cappellari et al. 2011), CALIFA (Walcher et al. 2014), SAMI (Croom et al. 2012), and MaNGA (Bundy et al. 2015), have provided kinematic maps as well as stellar age and metallicity maps for thousands of nearby galaxies. By applying an orbit-superposition dynamical model to the IFU data (Cappellari et al. 2007; van den Bosch et al. 2008; Zhu et al. 2018a), it has been possible to infer the internal orbital structure of individual galaxies within the data coverage. Two parameters are usually used to characterize the orbits in these models: the circularity, λ_z, which represents the angular momentum of the stars, and the time-averaged radius, r, as a proxy for binding energy. The distribution of stellar orbits in a galaxy can thus be represented by the probability density of orbits in the phase-space of circularity, λ_z, versus radius, r. The multiple structures of galaxies can also be identified within these phase-space planes with, for example, the subdivision of the stellar component in a galaxy into cold (λ_z > 0.8), warm (0.2 < λ_z < 0.8), hot (|λ_z|< 0.2), and counter-rotating (λ_z < −0.2) orbits.

In recent years, the stellar orbit circularity distribution within 1 R_e has been obtained for 300 CALIFA galaxies across the Hubble sequence (Zhu et al. 2018b) and for a few hundred early-type galaxies in MaNGA (Jin et al. 2020). These works show that the fraction of dynamically hot orbits increases from low-mass to high-mass galaxies, consistently with an increase in the mass fraction of the traditional bulge (Weinzirl et al. 2009). The IFU-based orbital decomposition of large numbers of galaxies is now achievable (Zhu et al. 2018c). However, the formation of these orbital components has not yet been quantitatively connected to galaxies’ merger histories.

The details of merger histories have been (partially) uncovered in the closest galaxies, where single stars in the stellar halos can be resolved, for example in the case of the Milky Way (Helmi et al. 2018; Belokurov et al. 2018; Helmi 2020), the Andromeda galaxy (D’Souza & Bell 2018a; D’Souza & Bell 2018b), and NGC 5128 (Rejkuba et al. 2011), or where individual globular clusters in the halos are detected (Forbes et al. 2016; Beasley et al. 2018; Kruijssen et al. 2019), or where streams have been detected (Merritt et al. 2016; Harmsen et al. 2017). Recently, studies have attempted to constrain the global ex situ fractions (Boecker et al. 2020; Davison et al. 2021) or the mass of satellite mergers (Pinna et al. 2019; Martig et al. 2021) with stellar population distributions obtained from IFU data in some edge-on lenticulars (Poci et al. 2019, 2021). However, so far the assembly and merger histories of most galaxies remain hidden, including in the nearby Universe.

From a numerical perspective, recent and current cosmological hydrodynamical simulation projects (Vogelsberger et al. 2020), such as IllustrisTNG¹ and EAGLE², have been able to reproduce large numbers of galaxies across the morphological spectrum with well-resolved structures (Correa et al. 2017; Du et al. 2019; Pillepich et al. 2019; Pulsoni et al. 2020). For example, the TNG100 and TNG300 simulations of the IllustrisTNG project produce galaxy morphological structures, such as galaxy concentration and bulge strength, that statistically match observational constraints for galaxies with M_* ≳ 10^9.5 M_⊙ at low redshift (Rodriguez-Gomez et al. 2019), greatly improving upon the analog outcome of the original Illustris simulation³. The stellar sizes of TNG100 and TNG300 galaxies have been shown to be in good agreement with observations at 10⁹ ≲ M_* ≲ 10¹¹ M_⊙, to within 0.1−0.2 dex across the full mass range (Genel et al. 2018). With the higher resolution of TNG50, galactic structures, in particular the stellar disk scale height and the stellar size of low-mass objects (Pillepich et al. 2019), match the observations even better (Zanisi et al. 2021). Similarly, the EAGLE simulation has been shown to successfully produce galactic structures, such as bulge sizes, that are similar to observations in galaxies with M_* ≳ 10¹⁰ M_⊙ (Lange et al. 2016); it has also produced galaxy stellar sizes that match observations well in galaxies with 10¹⁰ ≲ M_* ≲ 10¹¹ M_⊙ and are about 0.2 dex larger for lower or higher galaxy masses (Wang et al. 2020), at least when comparisons are done at face value (see de Graaff et al. 2022 for a discussion).

With the simulation data, and thanks to the realism of their simulated galactic structures, it has also been possible to trace back the formation of galactic structures in great detail (Lagos et al. 2018; Obreja et al. 2019; Trayford et al. 2019; Du et al. 2021; Pulsoni et al. 2021), including adopting dynamical decomposition methods similar to those applied to observational data. For example, by choosing the same definition of cold, warm, and hot components, Xu et al. (2019) systematically compared the orbital structures of TNG100-simulated galaxies to CALIFA galaxies and showed that the orbital fractions as a function of stellar mass in TNG100 are quantitatively consistent with the observed ones.

The availability of flexible dynamical decomposition methods that are applicable to large numbers of observed galaxies and the realism of current cosmological galaxy simulations means it is now possible to quantitatively combine the two to address the question of whether we can identify galactic structures in present-day galaxies that (1) are physically connected to their past merger events and (2) can be directly compared across simulation outputs and observations, given the unavoidable limitations on both ends. In this paper we specifically address this goal by focusing on the hot inner stellar halos of galaxies defined by stars with λ_z < 0.5 and 3.5 kpc < r < 2 R_e.

In the Milky Way, stars on highly radial orbits in the vicinity of the Sun have been associated with the Galaxy’s inner stellar halo and have been used as a probe of its ancient massive merger event(s) – such as the so-called Gaia Enceladus or Gaia Sausage progenitor (Belokurov et al. 2018; Helmi et al. 2018). Such structures in Milky Way-like galaxies can also arise in galaxy formation simulations (e.g., Grand et al. 2020). In this paper we address whether dynamical structures that are similar to the inner stellar halo of the Galaxy and that directly trace the ancient mergers generally exist across the mass and galaxy-type spectrum and whether they can be detected in external galaxies with current IFU data.

In practice, in this paper we analyze thousands of galaxies from the IllustrisTNG and EAGLE simulations, identify their hot inner stellar halo components in a way that is also consistent and achievable with current IFU observations, and – via the information provided by the merger trees of the simulations – look for quantitative connections between the inner stellar-halo mass and properties of the merger and stellar assembly histories. In particular, we focus on results from the TNG50 simulation and, after comparing them to those from other IllustrisTNG and EAGLE runs, use it to identify the main formation channels of the hot inner stellar halos of simulated galaxies. Our analysis is tailored toward the information that is obtainable from the Fornax3D survey (Sarzi et al. 2018), which uses MUSE/VLT and is designed to cover at least the inner 2 R_e of galaxies. In a companion paper, we will show that the stellar orbit distribution, as well as the stellar age and metallicity distributions, of each Fornax3D galaxy can be obtained through the population-orbit superposition method of Zhu et al. (2020), hence significantly improving upon the orbit characterization within 1 R_e with CALIFA or MaNGA observations and at last starting to probe into the stellar halos.

The paper is organized as follows: in Sect. 2 we describe the adopted cosmological simulations; in Sect. 3 we investigate the orbital decomposition of the simulated galaxies and define their hot inner stellar halos; in Sect. 4 we present the main results of the paper, which are then discussed and summarized in Sects. 5 and 6, respectively.

2. Cosmological galaxy simulations

2.1. TNG50 and IllustrisTNG

Among the existing cosmological hydrodynamical simulations for the formation and evolution of galaxies, the IllustrisTNG simulations (hereafter TNG; Springel et al. 2018; Marinacci et al. 2018; Naiman et al. 2018; Pillepich et al. 2018a; Nelson et al. 2018) have been particularly successful in reproducing a broad range of observational findings (Nelson et al. 2019a). These include the galaxy mass-size relation at 0 < z < 2 already mentioned in the Introduction (Genel et al. 2018), but also the gaseous and stellar disk sizes and heights (Pillepich et al. 2019), galaxy colors, the stellar age and metallicity trends at z ∼ 0 as a function of galaxy stellar mass in comparison to SDSS results (Nelson et al. 2018), and resolved star formation in star-forming galaxies (Nelson et al. 2021), as well as the characteristics of the stellar orbit distributions from the CALIFA survey (Xu et al. 2019) and the kinematics of early-type galaxies in comparison to data from ATLAS-3D, MaNGA, and SAMI (Pulsoni et al. 2020).

All the details of the numerical and physical model implemented in the TNG simulations are given in the method papers by Weinberger et al. (2017) and Pillepich et al. (2018b). Furthermore, the suite comprises three flagship runs: TNG50, TNG100, and TNG300, with a cosmological volume of (50 Mpc)³, (100 Mpc)³, and (300 Mpc)³ and a stellar particle resolution of 8.5 × 10⁴ M_⊙, 1.4 × 10⁶ M_⊙, and 1.1 × 10⁷ M_⊙, respectively (Nelson et al. 2019a). The basic information of these simulations is given in Table 1.

For most of the analysis in this paper, we use data from the TNG50 run (Pillepich et al. 2019; Nelson et al. 2019b). It is uniquely well suited for our study as (1) it reaches the numerical resolution typical of many zoomed-in simulation projects – with stellar particles of 8.5 × 10⁴ M_⊙, so that the inner structural details of galaxies, such as thin disks, thick disks, bulges, and halos, are fairly well resolved – and (2) it still offers a large sample of galaxies, including massive ones, in a representative cosmological volume that encompasses diverse environments.

In the TNG simulations, halos and subhalos, and thus the basic properties of galaxies, are identified using the friends-of-friends (FoF) and SUBFIND algorithms (Springel et al. 2001; Dolag et al. 2009). Based on SUBFIND, there are 548, 4002, and 50291 galaxies in the TNG50, TNG100, and TNG300 boxes at z = 0, respectively, with 10^10.3 < M_* < 10^11.6 M_⊙.

The histories of all galaxies in the simulations can be followed via the merger trees, as constructed by the SUBLINK-GAL code based on the baryonic component of subhalos (Rodriguez-Gomez et al. 2015). In this algorithm, each galaxy is assigned a unique descendant. We can identify for each galaxy its “mergers” – secondary galaxies with a well-defined “infall” time that merge with the primary (i.e., do not exist after “coalescence” as individually identified SUBFIND objects). For each galaxy at a given time we can also identify its “satellites” – galaxies with a well-defined infall time that orbit or fly by the primary prior to their possible ultimate merging.

Thanks to the merger trees, we can also distinguish between in situ and ex situ (i.e., accreted) stars (Pillepich et al. 2018b; Rodriguez-Gomez et al. 2019). For any galaxy at z = 0, we call its stars “in situ” if they formed in a progenitor that belongs to the main-progenitor branch of the galaxy, and “ex situ” otherwise. Ex situ stars are hence stellar particles stripped and consequently accreted into a galaxy: the ex situ stars of a galaxy are mostly accreted from its mergers (i.e., from objects that are already merged at the time of inspection); however, a fraction of ex situ stars can also arise from the stripping of orbiting and not-yet destroyed satellites.

2.2. EAGLE

The EAGLE (Schaye et al. 2015; Crain et al. 2015) simulations have also been shown to successfully reproduce a range of observations of galactic properties, including the galaxy stellar mass function, the Tully–Fisher relation, and the galaxy mass-size relation. As mentioned in the previous section, the galaxy sizes as a function of stellar mass generally agree with the SDSS results, with differences smaller than 0.2 dex at 10¹⁰ ≲ M_* ≲ 10¹¹ M_⊙. Galactic structures, such as disks and bulges, are well resolved, and the Hubble sequence is in place (Lange et al. 2016; Trayford et al. 2019).

In this work we use publicly available data (McAlpine et al. 2016) of the fiducial EAGLE simulation with a cosmological volume of about (100 Mpc)³, with a stellar particle mass of ∼10⁶ M_⊙ and a Plummer equivalent gravitational softening of 0.7 in proper kpc at redshift z = 0. The galaxy formation model and numerical code that produced the EAGLE outcome differ substantially from those of TNG. Here, we work with EAGLE halos and subhalos as identified by the same FOF group finder and SUBFIND algorithms adopted for the TNG simulations, based on which we also ran the same merger tree code as for TNG (Nelson et al. 2019a). The analysis of the EAGLE and TNG simulations is thus identical in methodology, including the merger and assembly definitions. There are 2049 galaxies in the EAGLE box with 10^10.3 < M_* < 10^11.6 M_⊙ at z = 0.

2.3. TNG and EAGLE galaxies adopted in this work

We use TNG50 for all the analysis throughout the paper. TNG100, TNG300, and EAGLE are only used to check if our results are robust against different galaxy formation and feedback models and across a range of numerical resolution.

We selected simulated galaxies based on their galaxy stellar mass and stellar size. The left panel of Fig. 1 shows the stellar size-mass plane for TNG50 at z = 0. We define galaxy stellar mass as the total mass of all stars that are gravitationally bound to a galaxy according to the SUBFIND algorithm. For TNG50 galaxies, we calculated the circularized projected half-light sizes in the z band, R_{e, z band}, to make them comparable with results from SDSS galaxies (Shen et al. 2003). Most galaxies are late-type at M_* ∼ 10¹⁰ M_⊙ and are early-type at M_* > 10¹¹ M_⊙. The sizes of TNG50 galaxies are generally consistent with the stellar size-mass relation of SDSS galaxies at M_* > 10¹⁰ M_⊙, as already quantified in depth by, for example, Zanisi et al. (2021).

Fig. 1. Stellar mass and size of the simulated galaxies analyzed in this paper. Left and central panels: total stellar mass, M*, versus projected half-light size in the z band, Re, z band, of TNG50 galaxies at z = 0. The dashed black curve is the median size of TNG50 galaxies as a function of M*. The red and blue curves show the median sizes of early-type and late-type SDSS galaxies (Shen et al. 2003) for reference. TNG50 galaxies with 1010.3 < M* < 1011.6 M⊙ are color coded by their stellar ex situ fraction (see also the zoomed-in version in the middle panel). The ex situ fraction systematically increases with increasing stellar mass and size. However, there are still large variations among galaxies of similar mass and size. We further divide the TNG50 samples into nine subsamples as indicated by the dashed squares labeled g1–g9 from left to right and from top to bottom. Right panel: median stellar sizes as a function of M* for TNG50, TNG100, TNG300, EAGLE, and Illustris galaxies at z = 0. Here we use the 3D half-mass size, r1/2,  3D,  mass, for all galaxies. The vertical dashed lines mark the mass regions we focus on in this paper.

The TNG50 galaxies of Fig. 1 are also color coded by their stellar ex situ fraction (f_ExSitu = M_*, ExSitu/M_*) following the definition and methods used earlier (Rodriguez-Gomez et al. 2016; Pillepich et al. 2018a). The stellar ex situ fraction of galaxies increases with increasing stellar mass and size, but still with a large scatter in galaxies of similar mass and size. This result confirms a similar trend predicted by the EAGLE model (Davison et al. 2020).

Throughout the paper, we focus on galaxies with stellar masses in the range 10^10.3 − 10^11.6 M_⊙, irrespective of whether they are centrals or satellites. We further divide the TNG50 samples into nine subsamples. We first made three bins of stellar mass, then three bins in size for each mass bin. The sizes are divided as $R_{\mathrm{e}}>\overline{R_{\mathrm{e}}}+\sigma (R_{\mathrm{e}})/2$, $\overline{R_{\mathrm{e}}}-\sigma (R_{\mathrm{e}})/2<R_{\mathrm{e}}<\overline{R_{\mathrm{e}}}+\sigma (R_{\mathrm{e}})/2$, and $R_{\mathrm{e}}<\overline{R_{\mathrm{e}}}-\sigma (R_{\mathrm{e}})/2$, and thus the numbers of galaxies in the three size bins are similar. Here $\overline{R_{\mathrm{e}}}$ is the median and σ(R_e) is the 1σ scatter of R_e in each mass bin. The nine subsamples are labeled as g1−g9 from left to right and from top to bottom. In the right panel of Fig. 1, we compare the galaxy size-mass relations for TNG50, TNG100, TNG300, and EAGLE. Here we simply use the 3D half-mass radius, r_{1/2,  3D,  mass}. TNG50, TNG100, TNG300 and EAGLE galaxies have similar sizes at M_* < 10¹¹M_* and match the observations well, whereas TNG100, TNG300, and EAGLE galaxies have sizes about 0.2 dex larger than observations at the high-mass end (Genel et al. 2018).

The sizes of Illustris galaxies (Vogelsberger et al. 2014a,b; Genel et al. 2015; Nelson et al. 2015) are also shown for comparison (red curve). As noted in previous papers (Pillepich et al. 2018b; Rodriguez-Gomez et al. 2019), Illustris galaxies have larger sizes at M_* < 10¹¹M_* in comparison to the other simulations and to observations: in fact, an improvement on galaxy sizes was one of the goals of the TNG project (Pillepich et al. 2018b). As realistic galactic structures are a key requirement for us to define the hot inner stellar halo in a coherent way (Rodriguez-Gomez et al. 2019), we do not not use the original Illustris simulation to extract the main results of this paper.

3. Dynamical decomposition of galaxies

We aim to separate the stellar constituents of the simulated galaxies into kinematically motivated “morphological” components, such as disk, bulge, and hot inner stellar halo, to then focus on the formation channel of the hot inner stellar halo and its information content.

3.1. Orbital properties

We characterize the stellar orbits of simulated galaxies with two parameters, the time-averaged radius, r, and circularity, λ_z. We measured these values with two complementary approaches.

In the first method, the radius and circularity of each stellar particle were calculated as the averages of the corresponding quantities for particles stored with equal time steps along the orbits. For the simulated galaxies, we have the instantaneous values of the stellar particles at z = 0, but we can integrate the orbits within the final potential. Namely, we froze the stellar particles of a simulation at z = 0, constructed the simulated potential by adding up stars, dark matter, and gas mass distribution, and integrated their orbits in the potential based on their instantaneous positions and velocities. We then calculated the time-averaged λ_z and r along the orbits and took these as the particle’s orbital properties. By doing this, all the particles on box orbits, which could instantaneously have a large variety of λ_z values, will have λ_z ≃ 0 as the time-averaged values along their orbits.

In a second approach, we adopted an approximate method (Zhu et al. 2018b), whereby averaging is done in phase-space. In practice, we assumed stellar particles with similar energy, E, and angular momentum, L_z, to be on similar orbits: we hence measured the average r and λ_z across such particles from the simulation and adopted them as the orbital r and λ_z values. This method cannot return λ_z = 0 for all particles on box orbits, but it does narrow their λ_z distribution.

We integrated the stellar orbits (i.e., followed the first method) using the code of AGAMA (Vasiliev 2019) for four TNG50 case-study galaxies whose evolution and formation are studied in great detail (see Sect. 5.2), so throughout the paper their stellar orbit properties are time-averaged values. The integration takes ∼1 CPU hour for each galaxy. For the statistical study of the large galaxy samples from the TNG and EAGLE simulations, we adopted the second method, which is good enough for our four-component separation (see Fig. A.1). For the four case-study galaxies from TNG50, we can in fact obtain galactic components via both methods, and we find consistent results. Throughout the paper, the radius of the orbit, r, refers to the time-averaged radius or the equivalent phase-space-averaged radius.

3.2. Separation into four stellar components

3.2.1. General separation

In Fig. 2 we show the stellar orbit distribution, p(r, λ_z), the stellar age distribution, t(r, λ_z), and the stellar metallicity distribution, Z(r, λ_z), in the phase-space of circularity, λ_z, versus radius, r, of a simulated galaxy from TNG50 at z = 0 (Subfind ID 468590). We use the distributions in Fig. 2 to visualize the orbital decomposition that we use throughout the paper.

Fig. 2. Definition of hot inner stellar halo advocated for and adopted in this paper. We show the stellar orbit distribution, p(r, λz), stellar age distribution, t(r, λz), and metallicity distribution, Z(r, λz), in the phase-space of circularity, λz, versus radius, r, of a simulated galaxy, TNG50 468590. We define as a hot inner stellar halo the galactic component whose stars have orbits with rcut < r < rmax and λz < λz, cut, where rcut, rmax, and λz, cut are, respectively, 3.5 kpc, max(7 kpc, 2 Re), and 0.5.

We decomposed a galaxy into four components based on the structures in the phase-space of λ_z versus r. We first notice that, in the case of “disky” galaxies, dynamically cold orbits with λ_z > 0.8 are usually younger and more metal rich than the other orbits. In the simulated galaxies, these orbits represent a rather well-defined thin stellar disk.

For the remaining non-disk orbits, we identified the following cuts to separate hot inner stellar halos from possible bulges and “thick” disks. Namely, we defined the hot inner stellar halo by selecting the orbits with r_cut < r < r_max and λ_z < λ_z, cut; throughout this paper, r_cut is equal to 3.5 kpc, r_max is chosen based on the maximum radius covered by IFU observations in real galaxies, and λ_z, cut = 0.5.

In practice, choosing r_cut means identifying a way to separate between the bulge and the hot inner stellar halo. We discuss this choice in detail in the next subsection.

Due to the limited observational data coverage, we can only obtain the stellar orbit distribution of real galaxies out to certain galactocentric distances. In the Fornax3D survey (Sarzi et al. 2018), the data coverage extends to ∼2 R_e for most galaxies. Here we chose r_max = 2 R_e for galaxies with R_e ≥ 3.5 kpc, whereas for galaxies with R_e < 3.5 kpc we chose r_max = 7 kpc: these choices imply that we can focus on what people may identify or call the “inner” stellar halo of a galaxy, rather than the halo component that extends out to a distance of ∼100 kpc.

We further distinguished between hot and warm orbits beyond the bulge radial cut by imposing a separation at λ_z, cut = 0.5. The idea was to distinguish between halo stars and stars in warmer orbits that could be identified as thick disk stars and hence possibly exhibit more complex and mixed origins. In the discussion, we quantitatively show that the choice of λ_z, cut for the identification of hot orbits is not important so long as λ_z, cut ≲ 0.6 − 0.7.

In summary, we separated each galaxy into four components: (i) the disk (λ_z > 0.8 and r < r_max), (ii) the bulge (λ_z < 0.8 and r < r_cut), the warm component (0.5 < λ_z < 0.8 and r_cut < r < r_max), and the hot inner stellar halo (λ_z < 0.5 and r_cut < r < r_max). Here, r_cut = 3.5 kpc and r_max = 2 R_e for galaxies with R_e > 3.5 kpc, and r_max is fixed at 7 kpc otherwise.

3.2.2. Choice of r_cut

In simulated and observed galaxies, we find that non-disk stellar orbits are high density and metal rich in the inner regions and become progressively lower density and more metal poor toward larger galactocentric distances. However, how to separate between the two components, the bulge and the hot inner stellar halo, and thus how to choose r_cut, is not immediately obvious. In the following, we provide arguments in support for our choice, r_cut = 3.5 kpc.

In the top panels of Fig. 3 we show the probability density distributions in panel a and the average metallicities in panel b as a function of radius, r, of the dynamically hot orbits (i.e., with λ_z < 0.5) in a selection of TNG50 galaxies (one galaxy, one curve). In particular, the diagnostics here are shown for galaxies in the group g6 (see Fig. 1). As can be seen, the probability densities, p_hot, first decrease sharply from a peak at small radii (< 5 kpc) and then flatten. On the other hand, the metallicities, Z_hot, have their maximum at the very center of galaxies, decrease sharply at first, and then flatten. Each galaxy curve is color coded by the ex situ stellar mass fraction, and the probability density distributions (p_hot) are normalized such that the total probability of all stars within 2 R_e equals unity. At a similar galaxy mass, galaxies with a higher stellar ex situ fraction have a higher fraction of hot orbits (lower fraction of disk), whereas the central density peaks are similar and galaxies with a higher ex situ fraction exhibit a higher plateau at r ≳ 3.5 kpc in hot orbit densities and at higher metallicity values.

Fig. 3. Properties of the stellar orbit distributions of simulated galaxies and physical justification of the definition of hot inner stellar halo. Top: probability density distribution in panel a and the average metallicity in panel b, as a function of radius, r, obtained by summing up all the hot orbits with λz < 0.5 from a selection of TNG50 galaxies at z = 0, namely those in group g6 of Fig. 1. Each galaxy’s curve is color coded by the galaxy’s stellar ex situ fraction, as indicated by the colorbar. The vertical line indicates the radius rcut = 3.5 kpc, where both the probability density and metallicity profiles flatten. Middle: distributions of the radii where the probability density distribution of hot orbits peaks, r|max(phot), for galaxies with M* > 1010.3 M⊙ in the TNG50, TNG100, TNG300, EAGLE, and Illustris simulations. Most galaxies in TNG50 have r|max(phot) < 1.5 kpc, consistent with real galaxies from the Fornax3D survey. Galaxies in TNG100, TNG300, and EAGLE have slightly larger bulge sizes but still with r|max(phot) < 3.5 kpc. Bottom: panel d: cumulative distribution function of the hot orbits with λz < 0.5, CDFhot as a function of radius, r, for a selection of TNG50 galaxies. Symbols and selection are the same as in the top panels. The green curve represents the control model constructed by the galaxies with the lowest stellar ex situ fraction. Bottom: panel e: deviation of CDFhot from the control model, CDFhot, control. The vertical line indicates the radius rcut = 3.5 kpc, which we adopt throughout as the inner boundary of the hot inner stellar halo.

We extended this characterization to all galaxies with 10^10.3 < M_* < 10^11.6 M_⊙ in the TNG50, TNG100, TNG300, and EAGLE simulations by measuring the peaking radius, r|_max(phot), for each galaxy (i.e., the radius where p_hot reaches the maximum). Again for comparison, we also show results from Illustris galaxies. In panel c of Fig. 3 we can see that most TNG50 massive galaxies have r|_max(phot) < 1.5 kpc: this is consistent with measurements of real galaxies of similar mass and size from the Fornax3D survey, as we will show in the companion paper. On the other hand, galaxies in other simulations are less dense in the very inner regions, which is shown by the larger r|_max(phot) values in Fig. 3: this indicates somewhat overly large bulge sizes, which can be due to a combination of a larger representation of more massive galaxies. This, in turn, is due to the larger volumes of the simulations in comparison to TNG50, lower numerical resolutions, larger force softening length, and/or different galaxy model prescriptions. In particular, a large fraction of Illustris galaxies have r|_max(phot) > 3.5 kpc: this is not consistent with real galaxies (see the companion paper for galaxies from the Fornax3D survey), and therefore we do not include Illustris galaxies in any further analysis. In contrast, TNG100 galaxies have r|_max(phot) similar to TNG50 but shifted slightly to larger values, as would be expected given the lower resolution.

To quantitatively determine a physically motivated transition radius, r_cut, we show in panel d of Fig. 3 the cumulative distribution function of the hot orbits, CDF_hot, of the TNG50 galaxies of the top panels. We also construct a control model by using the galaxies with the lowest stellar ex situ fraction (f_ExSitu − min(f_ExSitu) < 0.05) in the subsample. In panel e we show the deviation of the cumulative distribution function of each galaxy from the control model, CDF_{hot, control} − CDF_hot, again color coded by the stellar accretion fraction of the galaxy.

The underlying idea of this comparison is as follows: we postulate that the dynamically hot orbits of a galaxy would form a pure bulge if the galaxy had no mergers, whereas we call the component that would be formed by mergers the hot inner stellar halo. The “pure” bulge model is unknown in detail, but we approximate it by the control model that is constructed by galaxies with the smallest ex situ fraction in the subsample. Therefore, we take the deviation from the control model to indicate a sort of perturbation caused by the mergers (i.e., the contribution to a hot inner stellar halo component).

We hence defined the transition radius between the bulge and the hot inner stellar halo as the radius (r|_{max(CDF deviation)}) where the deviation, CDF_{hot, control} − CDF_hot, reaches its maximum. This is distributed closely around 3.5 kpc for TNG50 galaxies in groups g2, g5, g6, and g9, though with larger scatter in g3, and it is slightly smaller and with larger scatter in g1, g4, g7, and g8 (see Fig. B.1). Galaxies in each group have a range of sizes, R_e, while r|_{max(CDF deviation)} is distributed closely around a physical size of ∼3.5 kpc. The bulges in these galaxies have similar physical sizes, consistent with the results from Du et al. (2020); the bulge size does not scale with R_e, which is instead mostly set by the buildup of the outer disk or the halo.

To keep the analysis simple and easy to compare with real observations, we chose r_cut = 3.5 kpc universally across all galaxies. The choice of r_cut = 3.5 is not perfect for the galaxies from the other simulations, as they exhibit larger bulges and thus somewhat larger transition radii. However, a bulge to halo separation at 3.5 kpc still applies as a well-motivated choice for TNG100, TNG300, and EAGLE galaxies with 10^10.3 < M_* < 10^11.6 M_⊙, as for most of them the maximum of p_hot is reached at r|_max(phot) < 3.5 kpc, and thus the hot-orbit distributions p_hot are declining at 3.5 kpc.

4. Information content of the hot inner stellar halo

We took all the galaxies with 10^10.3 < M_* < 10^11.6 M_⊙ in the TNG50, TNG100, TNG300, and EAGLE simulations at z = 0. For each of them, we measured the mass of their hot inner stellar halo, M_{*, halo(r < 2 Re)}, by adopting the orbital decomposition described in Sect. 3.2: we did so irrespective of their global morphology, kinematics, star formation properties, and environments.

We also kept track of the merger and stellar accretion histories of each simulated galaxy and summarized them by means of statistics, such as the total ex situ stellar mass. This is the total amount of stellar mass accreted via mergers and via the stripping of satellites (see Sect. 2.1) and the galaxy stellar mass of the most massive progenitor ever merged with.

In the following we show that the mass of the hot inner stellar halo alone is sufficient to infer selected summary statistics of the merger and assembly history of individual galaxies.

4.1. Mass fraction of the hot inner stellar halo

In Fig. 4 we show how the mass fraction of the hot inner stellar halo is connected to the galaxy’s stellar ex situ fraction, f_ExSitu. We define the mass fraction of the hot inner stellar halo compared to the total stellar mass within the maximum radius of r_max = 2 R_e as f_halo = M_{*, halo(< 2 Re)}/M_{*(< 2 Re)}.

Fig. 4. Connection between the hot inner stellar halo and the ex situ mass fraction according to TNG50. We show the mass fraction of the hot inner stellar halo compared to the total stellar mass within the maximum radius of rmax = 2 Re (fhalo = M*, halo(< 2 Re)/M*(< 2 Re)) versus the stellar ex situ fraction (fExSitu = M*, ExSitu/M*) for galaxies with M* > 1010.3 of TNG50 at z = 0. Here we exclude 37 objects identified with ongoing mergers or those with counter-rotating thin disks. Galaxies with size Re below, around, and above the average are indicated by dots, asterisks, and crosses, respectively: they are color coded based on stellar mass, with increasing stellar mass from black, to blue, to red. Galaxies with larger stellar masses and sizes exhibit systematically larger ex situ fractions as well as larger mass fractions in the hot inner stellar halo. The Pearson correlation coefficient ℛ(fhalo, fExSitu) = 0.86 of the whole sample is annotated. The black line is a linear fit to the data points, while the 1σ scatter of the data points against the linear fit is σ = 0.07.

By checking the stellar images of all TNG50 galaxies in the mass range shown in the figure (with M_* = 10^{10.3 − 11.6} M_⊙), we identified 37 that are undergoing a merger or exhibit a counter-rotating thin disk at the time of inspection (z = 0). Galaxies with ongoing mergers are not in dynamical equilibrium, making a meaningful definition of the hot inner stellar halo impossible. On the other hand, counter-rotating thin disks are usually formed in situ from misaligned gaseous disks, regardless of mergers. These 37 galaxies significantly deviate from the correlation and are excluded from the figure.

Mass fraction of the hot inner stellar halo, f_halo, is strongly correlated with the galaxy’s stellar ex situ fraction, f_ExSitu. We get a Pearson correlation coefficient of ℛ(f_halo, f_ExSitu) = 0.86 with the data points in the figure, and it decreases to 0.82 when the 37 objects are included. The black line is a linear fit to the data points, and we calculate that the 1σ scatter of the data points against the linear fit is σ = 0.07.

In Fig. 4 we show that galaxies with larger stellar mass, M_*, and size, R_e, have systematically larger f_ExSitu and f_halo. Thus, the variation in f_ExSitu across the mass-size plane (as shown in Fig. 1) is encoded in the variation in f_halo. Moreover, galaxies with similar stellar mass and size still span a wide range in f_ExSitu and f_halo. The relation we show here reflects the physical connection between the increase in the stellar ex situ fraction and the growth of the hot inner stellar halo, regardless of the galaxy stellar mass and size.

4.2. Correlations between the hot inner stellar halo mass and the ancient accreted stellar mass

Figure 5 gives the main results of this paper. There, we show the total ex situ stellar mass, M_*, ExSitu (left), and stellar mass accreted from the most-massive galaxy ever merged by TNG50 galaxies, M_*, Ex1 (right), as a function of their hot inner stellar halo mass, M_{*, halo(r < 2 Re)}. Galaxies are color coded by their total stellar ex situ fraction⁴.

Fig. 5. Correlations between the mass of the hot inner stellar halo, M*, halo(r < 2 Re), and the total accreted stellar masses (left) and the stellar mass accreted from the most massive merger (right), according to TNG50. Left panel: M*, ExSitu includes all the ex situ stellar mass ever accreted from all past mergers and from the stripping of ongoing mergers, satellites, and flybys. Right panel: M*, ExSitu is the stellar mass accreted from the most massive merger. Here all galaxies with M* > 1010.3 from TNG50 at z = 0 are included, color coded by their stellar ex situ fraction. The thick black curve is the running median, and the thin dashed curves represent the ±1σ scatter. Black plus symbols mark galaxies with ongoing merging or those with counter-rotating disks. The Pearson correlation coefficients, ℛ, of the correlations are annotated in each panel. The dashed magenta line refers to 1:1. Some satellite mergers will deposit their stellar mass across a wide radius, with only a fraction of it within 2 Re, and so can lie above that line. On the other hand, the x axis shows a combination of in situ and ex situ stars such that the total mass can greatly exceed the total mass of the satellites in some galaxies with low ex situ fractions, causing them to lie below the line.

It is manifest from Fig. 5 that the hot inner stellar halo mass is tightly correlated with the total ex situ stellar mass, with Pearson correlation coefficient ℛ(M_{*, halo(r < 2 Re)}, M_*, ExSitu) = 0.91, as well as with the stellar ex situ fraction. The correlation with the stellar mass accreted from the most massive progenitor a galaxy has ever merged with is slightly weaker, with ℛ(M_{*, halo(r < 2 Re)}, M_*, Ex1) = 0.85.

The 1σ scatter of the correlations against the running medium is ∼0.5 dex at M_{*, halo(r < 2 Re)} < 10¹⁰ M_⊙. It becomes as small as ∼0.1 dex for M_*, ExSitu and ∼0.15 dex for M_*, Ex1 at the high-mass end explored here, whereby M_{*, halo(r < 2 Re)} ∼ 10¹¹ M_⊙ (i.e., a total galaxy stellar mass of 10^11.6 M_⊙).

The 37 galaxies that are undergoing a merger or exhibit a counter-rotating disk at the time of inspection (z = 0) are included in Fig. 5: they are typically outliers, most of them lying below the median trend. The Pearson correlation coefficients ℛ(M_{*, halo(r < 2 Re)}, M_*, ExSitu) and ℛ(M_{*, halo(r < 2 Re)}, M_*, Ex1) of the whole TNG50 sample increase to 0.92 and 0.86, respectively, when the 37 galaxies are excluded. Apart from this, most outliers lying below the correlations are galaxies with a stellar ex situ fraction f_ExSitu < 0.1.

It is important to highlight that the correlation between the hot inner stellar halo mass and the total ex situ stellar mass (or the mass of the most massive merger) is stronger than the one between the latter and the total stellar mass of a galaxy: this is the case even though all quantities of Fig. 5 correlate with a galaxy’s mass. This can be learned from Fig. 4, whereby it is manifest that the mass fraction of the hot inner stellar halo is highly correlated with the stellar ex situ fraction. We also quantify this in Appendix C: for comparison, for the same TNG50 galaxies, the Pearson correlation coefficient between a galaxy stellar mass and M_*, Ex1 reads 0.75 and the 1σ scatter varies from 0.6 to 0.3 from the low-mass to the high-mass end (see Fig. C.1).

In Fig. 6 we demonstrate that the existence of tight correlations is a prediction not exclusive to TNG50: we see that there are similar correlations for TNG50 (repeated from Fig. 5), TNG100, TNG300, and EAGLE galaxies. In the left panel, we show the relationships between the hot inner stellar halo mass of a galaxy and the total ex situ stellar mass, and on the right we give the relationship between the former and the accreted stellar mass from its most massive merger.

Fig. 6. Information content of the hot inner stellar halo across galaxy formation simulations. We show the correlations of M*, halo(r < 2 Re) with the total ex situ stellar mass, M*, ExSitu (left), and with the stellar mass accreted from the most massive merger, M*, Ex1 (right), for different galaxy simulations. The thick black, blue, yellow, and green lines are the median curves from TNG50, TNG100, TNG300, and EAGLE galaxies at z = 0, respectively, in the M* = 1010.3 − 11.6 M⊙ range; the thin curves are the ±1σ scatters. The Pearson correlation coefficients, R, from the four sets of simulations are labeled in the corresponding colors.

From these comparisons we can tell that results from models with different galaxy-formation physics and resolution agree to a good degree, particularly toward the high-mass end. In particular, the correlations from TNG100 galaxies match those from TNG50 perfectly, with very similar median and scatter throughout the mass range explored here. With lower resolution than TNG50, TNG100 galaxies have less-dense centers and larger bulge sizes than TNG50 galaxies, as discussed with Fig. 3. Here we use the same r_cut = 3.5 kpc for all galaxies to separate the hot inner stellar halo from the bulge, even though the actual transition radius of TNG100 galaxies may be somewhat larger. The fact that the correlation also holds for TNG100 galaxies is an indication that the precise choice for the inner boundary of the hot inner stellar halo is not crucial so long as the structural properties of the simulated galaxies are reasonably realistic.

At even lower resolution and with larger softening length, TNG300 galaxies have even larger bulge sizes than TNG100 galaxies (again, Fig. 3, middle panel). Our fixed choice of r_cut = 3.5 kpc causes an overestimate of the hot inner stellar halo mass, especially in galaxies with low ex situ fractions. The trends from TNG300 galaxies are consistent with those from TNG50 and TNG100 galaxies at M_{*, halo(r < 2 Re)} ≳ 10^10.2 M_⊙ ( i.e., galaxy stellar mass ≳10¹¹ M_⊙) but are steeper at lower masses.

EAGLE galaxies are the result of different feedback physics but have numerical resolution comparable to TNG100 galaxies and exhibit bulge sizes between those of TNG100 and TNG300 galaxies. The EAGLE correlations between the hot inner stellar halo mass and the ex situ stellar mass (total or from the most massive merger) are perfectly consistent with those from the TNG50 and TNG100 simulations at M_{*, halo(r < 2 Re)} ≳ 10^9.8 M_⊙ (i.e., galaxy stellar mass ≳10^10.6 M_⊙) but become steeper at M_{*, halo(r < 2 Re)} ≲ 10^9.8 M_⊙. Yet, the EAGLE predictions are within the 1σ scatters of the TNG50 and TNG100 predictions.

In fact, and importantly, all models predict tight correlations between the mass of the hot inner stellar halo and the total ex situ mass or the stellar mass of the most massive progenitor, with Pearson correlation coefficients larger than 0.86 and 0.8, respectively, across all simulation runs. The trends predicted by TNG50, TNG100, and EAGLE are consistent within the ±1σ scatter across the explored mass range.

4.3. Correlations for galaxies of different mass and size

In Fig. 1 we show that the ex situ stellar mass fraction of galaxies increases with increasing stellar mass and stellar size. In fact, even when the galaxy population is divided into subsamples at fixed galaxy stellar mass and varying size (see the nine groupings of Fig. 1), non-negligible variations in ex situ fraction, and hence merger histories, are manifest. Moreover, Fig. 5 shows that the relationship between the hot inner stellar halo mass and the total ex situ stellar mass is tighter the larger the galaxy mass. In the following, we quantify how the correlations uncovered in this paper change as a function of galaxy mass and stellar size.

In Fig. 7 we show the correlations of the hot inner stellar halo, M_{*, halo(r < 2 Re)}, versus M_*, ExSitu for the subsamples of TNG50 galaxies identified in Fig. 1: from less to more massive galaxies from left to right, from larger to smaller galaxies from top to bottom, with the galaxies in the middle row exhibiting median stellar sizes given their stellar mass. We performed the linear fitting and calculation of the Pearson correlation coefficient many times by bootstrapping. The medium and ±1σ uncertainty of the Pearson correlation coefficient, ℛ(M_{*, halo(r < 2 Re)}, M_*, ExSitu), and slope, p, of the linear fitting are labeled in each panel. The typical 1σ scatter of data points against the linear fit are also labeled.

Fig. 7. Correlation between M*, halo(r < 2 Re) and the total ex situ stellar mass, M*, ExSitu, as a function of galaxy stellar mass and stellar size. Here we show the relationships for TNG50 galaxies divided into nine groups by mass and size, as shown in Fig. 1, one grouping per panel, with galaxy stellar mass increasing from left to right and galaxy stellar size increasing from bottom to top. In the middle row, galaxies of median size for their mass are depicted. The black dots denote TNG50 galaxies without ongoing mergers, whereas the gray plus symbols indicate ongoing mergers. The thick black line is a linear fit to the black dots, and the thin black lines represent the ±1σ scatter. The median and ±1σ scatter of the whole sample as shown in Fig. 5 are reproduced as gray curves for reference. The medium and ±1σ uncertainty of the Pearson correlation coefficient, ℛ(M*, halo(r < 2 Re), M*, ExSitu), slope of the linear fit, p, and the typical 1σ scatter of the data points against the linear fit are labeled for each group. The four galaxies that are discussed as case studies in Sect. 5 are marked with colored symbols (in panels g6 and g9).

Figure 7 demonstrates that the information content of the hot inner stellar halo is maximal for the most massive and largest galaxies (i.e., the correlations are strongest in the subsamples with the most massive galaxies and/or those with a large size). The 1σ scatter is ∼0.12 dex in g2, g3, and g6: these values are similar to the scatter inferred from the whole galaxy sample at the high-mass end. On the other hand, the correlations are weaker in galaxies with lower masses and smaller sizes. The 1σ scatter increases to ∼0.25 dex in g1 and g5, and to 0.3−0.4 dex in g4, g7, g8, and g9 (i.e., mostly for galaxies below the average mass-size relation).

Figure 8 shows the analog correlation of the hot inner stellar halo mass, M_{*, halo(r < 2 Re)}, versus the stellar mass accreted from the most massive merger, M_*, Ex1, for each subsample of TNG50 galaxies. The correlations in most groups are slightly weaker when only the most massive merger is considered (Fig. 8 versus Fig. 7). The 1σ scatter is ∼0.2 dex in g2, g3, and g6 and increases to ∼0.4−0.5 dex in groups of galaxies with lower masses and smaller sizes.

Fig. 8. Correlation between the hot inner stellar halo mass, M*, halo(r < 2 Re), and the stellar mass accreted from the most massive merger, M*, Ex1, as a function of galaxy mass and size. Annotations are as in Fig. 7, and results are from TNG50.

Importantly, Figs. 7 and 8 show that the average relationships between the mass of the hot inner stellar halo and the two studied summary statistics of the merger history of galaxies hold practically unchanged across the investigated galaxy population, barring the 10^{10.3 − 10.6} M_⊙ galaxies with below-average stellar sizes. The slopes of the linear fits of the nine groups are consistent with one another within 1σ uncertainty. On the other hand, the constraining power of the hot inner stellar halo mass is highest for more massive and more extended galaxies.

There are still positive, albeit very weak, correlations between the total stellar mass, M_*, on the one side and M_*, ExSitu or M_*, Ex1 on the other, including in narrow bins of galaxy stellar mass and size. The correlations in M_{*, halo(r < 2 Re)} versus M_*, ExSitu (or M_*, Ex1) are much stronger.

5. Discussion

We want to determine the physical origin of the relationships quantified and uncovered in Figs. 5–8 as well as the origin of the lingering galaxy-to-galaxy variations, which are also at fixed galaxy mass and size. To address these questions, we delve into the simulation data, follow the formation of the hot inner stellar halo of selected TNG50 galaxies, and check whether other properties can be used to provide additional information into the past assembly history of a galaxy. Before that, however, we comment on how the correlations presented in this paper depend on the adopted operational definition of a hot inner stellar halo.

5.1. Dependence on the definition of hot inner stellar halo

In Sect. 3.2 we chose a definition of hot inner stellar halo such that the circularity of its stars (or stellar orbits) is smaller than a given value, λ_z, cut = 0.5. Namely, in our operational definition of (inner) stellar halo, we select only for hot orbits and exclude possibly “warm” and “cold” ones – as a reminder, the disk component is composed of cold orbits with circularity λ_z > 0.8. Here we show how the correlations discovered above depend on such a choice.

In Fig. 9 we show the Pearson correlation coefficient, ℛ(M_{*, halo(r < 2 Re)}, M_*, Ex1), as a function of λ_z, cut for TNG50 galaxies, with the latter varying between about 0.2 and 1. We study results for both the whole TNG50 sample (see the selection in Sect. 2.3) and the TNG50 galaxies binned by galaxy stellar mass and size (see groupings from Fig. 1). The correlation is dramatically stronger with λ_z, cut ≲ 0.8 for most galaxy subsamples. So long as λ_z, cut ≲ 0.6 − 0.7, the correlations are similarly informative to within 10 − 20%, irrespective of galaxy properties. The operational definition of the hot inner stellar halo adopted throughout this analysis (with λ_z, cut = 0.5) is a favorable one for the scientific questions of interest but is not fine-tuned.

Fig. 9. Dependence of the correlation between the masses of the most massive merger and of the hot inner stellar halo of a galaxy on the operational definition of the latter. We show the Pearson correlation coefficient between M*, Ex1 and M*, halo(r < 2 Re) as a function of λz, cut, the maximum level of orbital “coldness” adopted to define inner-halo stars. Namely, in the plot we adopt different definitions of hot inner stellar halo by imposing λz < λz, cut. Dashed curves in different colors indicate galaxies in subsamples g1–g9, and the solid black curve indicates galaxies of the whole sample.

5.2. Formation of the hot inner stellar halo

We illustrate the formation of the hot inner stellar halos with four simulated galaxies from TNG50: namely, galaxies with Subhalo ID at z = 0 reading 2, 468590, 117253, and 375073. These represent the high-mass end of the range under scrutiny, with galaxy stellar masses of 1.9, 1.4, 1.3, and 1.9 × 10¹¹ M_⊙, respectively. However, TNG50 2 and 468590 are on the average mass-size relation, with stellar sizes, R_e, measuring 6.5 and 7.0 kpc, respectively – they hence belong to subsample g6, as annotated in Figs. 7 and 8. Galaxies TNG50 117253 and 375073 lie below the average mass-size relation for their mass (grouping g9), with sizes of 3.1 and 5.2 kpc, respectively. They have a variety of accretion histories, but all have a substantial hot inner stellar halo at z = 0, between a few 10¹⁰ and 10¹¹ M_⊙. In the following, we generally refer to “massive mergers” as to those whose secondary’s stellar mass is of the same magnitude as the most massive one ever accreted (1:10); the rest are considered to be minor mergers.

The formation history of TNG50 468590 is shown in panels a-1 and b-1 of Fig. 10. That of TNG50 375073 is shown in panels a-2 and b-2.

Fig. 10. Formation of the hot inner stellar halo in TNG50 468590 (a-1 and b-1; TNG50 468590 is the same galaxy as illustrated in Fig. 2) and TNG50 375073 (a-2 and b-2). Panel a-1: stellar mass assembly history of TNG50 468590. The black dots indicate the stellar mass evolution of the main progenitor as a function of redshift, z. The dashed colored lines denote the mass assembly of the secondary galaxies that merged into the main progenitors directly and at one point reached a maximum stellar mass larger than 109 M⊙; the corresponding colored star symbols indicate the times when the maximum mass was reached. TNG50 468590 experienced a massive merger at z ∼ 1 with a stellar mass ratio of ∼1:1. The solid and dashed red vertical lines indicate the beginning and ending of this merger. Top row of panel b-1: morphology of the main progenitor before the major merger in the leftmost panel and the morphology at redshift z = 0 in the rightmost panel; in between we have the contribution to the z = 0 morphology from different origins with (from left to right) the stars from the main progenitor, accreted from the most massive merger, formed during the major merger, formed after the major merger ends, and accreted from subsequent minor mergers. We use co-moving distances in units of kpc (labeled as ckpc) for galaxies at z > 0. The surface density is in units of stellar mass per unit area kpc2. The stellar mass is normalized such that the total stellar mass of the galaxy at z = 0 equals unity within the figure coverage. The bottom row shows the orbital distribution of the corresponding stars in the phase-space of circularity, λz, versus radius, r, as derived from the simulation particle information. The probability density is in units of stellar mass per unit area in the phase-space. The dashed lines indicate our orbital-based division into four components as adopted for NGC 1380: disk, warm component, bulge, and hot inner stellar halo. It has a hot inner stellar halo mass of 3.4 × 1010 M⊙, and 90% of it was produced by the most massive merger event. The contribution from subsequent minor mergers is negligible. Panels a-2 and b-2: are similar but for TNG50 375073. TNG50 375073 had two massive mergers, one at z ∼ 2 and one at z ∼ 1. Panel b-2: we trace the structure formation during the most massive merger, which happened at z ∼ 1. It has a hot inner stellar halo mass of 3.0 × 1010 M⊙, and 84% of it was produced by the most massive merger event.

The assembly history of TNG50 468590 (top panel of Fig. 10) shows that this galaxy’s past is dominated by a massive merger with a stellar mass ratio of about 1:1 at redshift z ∼ 1. This galaxy has accreted stellar mass for a total of 6 × 10¹⁰ M_⊙ (ex situ fraction of 43%), with its most massive merger contributing 94% of it. In panel b-1 we show the galactic structure of the progenitor of TNG50 468590 before and after this merger. Before the merger, the main progenitor was dominated by a cold disk and a compact bulge with no substantial hot inner stellar halo. During the merger, both the main progenitor and the satellite were destroyed, with some stars settling into the bulge and others forming part of the hot inner stellar halo. Both the main progenitor and the satellite were gas rich. As such, 4 × 10¹⁰ M_⊙ of stars (comparable to the mass of the main progenitor at the time) formed during the merger, most of which contributed to the bulge and ∼15% were distributed into the hot inner stellar halo.

After the merger, new stars formed on dynamically cold disk-like orbits. This disk structure persisted until z = 0 because the quiescent merger history after this event results in only little dynamical heating and a minor contribution of these ex situ stars to the hot inner stellar halo. The hot inner stellar halo mass of TNG50 468590 at z = 0 is 3.4 × 10¹⁰ M_⊙: importantly, only ∼10% of the inner-halo stars were already there in the main progenitor before the massive merger started; ∼30% were in the bulge and/or disk of the main progenitor and redistributed to the hot inner stellar halo during the merger; ∼40% were brought in by the satellite galaxy; and 20% formed during the merger and were immediately distributed into the hot inner stellar halo regions. In total, about 90% of the hot inner stellar halo of this simulated galaxy was the product of the massive merger event.

Similar figures for TNG50 375073 in panels a-2 and b-2 of Fig. 10 show that TNG50 375073 has undergone two massive mergers, one of which contributed 58% of the total accreted stellar mass (i.e., 5.1 × 10¹⁰ M_⊙). This galaxy has a hot inner stellar halo mass of 3.0 × 10¹⁰ M_⊙ at z = 0. We also show in this case the galactic structure during the most massive merger, which started at z ∼ 1.3 with a stellar mass ratio of ∼1:0.5. At z = 1.3, before this merger started, ∼14% of the hot inner stellar halo stars were already in the main progenitor, which might be the result of the less massive merger that happened earlier on. During the most massive merger, ∼38% of the stars of the z = 0 inner halo were redistributed into the halo regions from the bulge and disk of the main progenitor; ∼34% were accreted from the satellite galaxy; and ∼12% were formed during the merger event. In total, 84% of the hot inner stellar halo mass of TNG50 375073 was induced by the most massive merger event.

The other two galaxies are characterized by a smaller fraction of the hot inner stellar halo mass produced by the most massive merger event. They are shown in Fig. 11.

Fig. 11. Formation of the hot inner stellar halo in the simulated galaxies TNG50 2 (a-3 and b-3) and TNG50 117253 (a-4 and b-4). The symbols are the same as in Fig. 10. Panel a-3: TNG50 2 has accreted three massive satellites, which are denoted by Ex1, Ex2, and Ex3, with decreasing stellar mass. Here, Ex3 and Ex2 were accreted simultaneously at z ∼ 1 and marked as the first merger, and Ex3 was accreted at z ∼ 0.1 and marked as the second merger. Panel b-3: we trace the structure formation during the two merger events. During the first merger, stars accreted from both Ex2 and Ex3, as well as stars formed during the merger, contribute to the hot inner stellar halo. After the first merger and before the second merger, a large number of stars were formed in a cold disk. The disk was then heated, becoming warm, during the second merger, and stars accreted from Ex1 contribute significantly to the hot inner stellar halo. Stars are rarely formed during and after the second merger as there is little gas left. TNG50 2 has a hot inner stellar halo mass of 4.3 × 1010 M⊙, of which 45% is a product of the most massive merger. Panel a-4: TNG50 117253 has accreted three massive satellites, which are still denoted by Ex1, Ex2, and Ex3, with decreasing stellar mass. We focus on the last two accretion events, Ex1 and Ex2, at 1 < z < 1.8. Panel b-4: structure formation during the accretion of Ex1 and Ex2. TNG50 117253 has a hot inner stellar halo mass of 1.5 × 1010 M⊙: the fraction produced by the most massive merger has an upper limit of 60%.

The galaxy TNG50 2 (panels a-3 and b-3 of Fig. 11) accreted a total of 9 × 10¹⁰ M_⊙ (ex situ fraction of 47%), with the most massive merger contributing 50% of it. We split its merger history into two parts: 1.5 < z < 0.9 when the second and third massive satellites (Ex2 and Ex3) were accreted, which we call the first merger, and 0.5 < z < 0.1 when the accretion of the most massive merger (Ex1) happened, which we call the second merger. This galaxy has a hot inner stellar halo mass at z = 0 of 4.3 × 10¹⁰ M_⊙, and only 5% of it was already in place before the first merger.

During the first merger, stars from Ex2 and from Ex3 are accreted to make up 29% and 5%, respectively, of the final mass of the hot inner stellar halo, whereas stars formed during this merger event contribute 14% of the hot inner stellar halo. After the first merger and before the second merger started, a cold disk was reformed. During the second merger, the satellite galaxy Ex1 deposited its stars, making up 40% of the hot inner stellar halo. The cold disk was heated to be dynamically warm, with a small part of it heated to be dynamically hot; furthermore, a small amount of stars were formed during the merger and contribute to the inner stellar halo. The combination of these two amounts to ∼5%. In this case, in total only about 45% of the hot inner stellar halo mass is a product of the most massive merger event.

Finally, in the case of TNG50 117253 (panels a-4 and b-4 of Fig. 11), the total accreted stellar mass is 3.2 × 10¹⁰ M_⊙, 35% of which comes from its most massive merger. Three massive satellites were accreted in quick succession between z = 2.5 and z = 1. By considering a snapshot of this galaxy at z = 1.8, when the first merger ended, we can see that 12% of the hot inner stellar halo mass was already in place at that time. Then another two satellites interacted with the main progenitor simultaneously, and it is hard to distinguish the effects of these two satellites. During these merger events, we find that ∼23% of the hot inner stellar halo stars were accreted from the most massive merger (Ex1); ∼15% were accreted from the second massive satellite (Ex2); and ∼40% were formed and induced during the mergers. At the end, TNG50 117253 had a total hot inner stellar halo mass of 1.5 × 10¹⁰ M_⊙. The fraction produced by the most massive merger event is less than 60%.

The four cases described above give a snapshot of the complex assembly history of galaxies and, in particular, of the formation processes of their hot inner stellar halo. Yet, in all four cases, the stars accreted from the most massive merger(s) favor highly radial orbits with λ_z ∼ 0. This is consistent with massive satellites on radial orbits being the most likely to deposit stars in the inner regions (Boylan-Kolchin et al. 2008). Still, the hot inner stellar halo is not only a result of the stripping of ex situ stars from incoming satellites and mergers, however massive.

Our analysis of TNG50 galaxies reveals that the stars in the hot inner stellar halo of galaxies have three main origin channels: (a) they can be accreted from the satellite(s); (b) they can be kicked out from the bulge and/or disk of the main progenitor by a massive merger; or (c) they can be formed during a massive merger. This is consistent with the “dual” halo found in the Milky Way (Haywood et al. 2018; Grand et al. 2020) and previously theorized for Milky Way-like galaxies by a number of numerical simulations (e.g., Zolotov et al. 2009; Pillepich et al. 2015).

A galaxy could undergo a few massive mergers during its history, and all in principle can contribute to the formation of its hot inner stellar halo. Based on the four example galaxies studied above, we can conclude that, depending on a galaxy’s merger and accretion history, its most massive merger event may directly (via ex situ stars) or indirectly (by inducing in situ star formation or in situ star displacement) contribute from 40% to 90% of the mass of the hot inner stellar halo.

5.3. Contribution of in situ and ex situ stars

From the study cases above, we have learned that the hot inner stellar halo as defined in Sect. 3.2 is a mixture of in situ and ex situ stars. Here we extend the analysis above to all galaxies in the TNG50 sample. In the top panel of Fig. 12 we show the probability density distribution of stars on dynamically hot orbits as a function of radius, separately for in situ (left panel) and ex situ (right panel) stars. Here we further divide the subsample of galaxies in g6 into five bins according to their ex situ fraction. Each curve represents the average p_hot (in situ or ex situ) of galaxies in each bin, color coded by its range of stellar ex situ fraction. The curves are normalized such that the total probability of all stars within 2 R_e equals unity. Thus, the galaxies with lower ex situ fractions have lower fractions of hot orbits, in situ and ex situ stars combined, due to higher fractions of disk orbits.

Fig. 12. Ex situ and in situ origin of the hot inner stellar halo of TNG50 galaxies. Top: probability density distribution of stars on hot orbits, for in situ (left) and ex situ (right) stars separately. Here we divide the subsample of TNG50 galaxies in g6 into five bins according to their ex situ fraction, fExSitu. Each curve represents the average phot (in situ or ex situ) of galaxies in each bin, color coded by its range of stellar ex situ fraction from 0−0.2 to 0.5−0.7, as indicated by the colored numbers. The curves are normalized such that the total probability of all stars within 2 Re equals unity. Bottom: stellar ex situ fraction of the galaxy, fExSitu, versus the ex situ fraction of stars in the hot inner stellar halo, fExSitu, halo (3.5 kpc < r < 2 Re and λz < 0.5). Each dot represents one galaxy in TNG50, color coded by galaxies stellar mass, M*, and crosses denote galaxies in subsample g6 (i.e., with a galaxy stellar mass in the range 1011 − 1011.6 and an average stellar size).

In galaxies with low ex situ fractions, the in situ hot stars are highly concentrated in the inner regions, with a low tail outside, whereas the ex situ stars are radially extended but with low density. In galaxies with higher ex situ fractions, the density tails of in situ stars at r > 3.5 kpc are higher. Combining this with what we learned from the example galaxies above, this means that the distributions of in situ stars on previous bulges and/or disks are altered by the mergers, with more massive mergers inducing more prominent tails of the in situ hot stars at a large radius. At the same time, ex situ stars also contribute more to the tail with higher ex situ fractions, although about half of the ex situ stars sink into the center and become part of the bulge. In galaxies with stellar ex situ fractions f_ExSitu ≳ 0.5, the density distributions of in situ and ex situ stars are similar, as they could have had major mergers with a mass ratio close to 1:1.

We note that a bulge with a similar size as that in galaxies with few ex situ stars formed after major mergers. The density of hot orbits peaks at a similar radius, with r|_max(phot) ∼ 1 kpc for galaxies with different ex situ fractions.

In the bottom panel of Fig. 12 we further quantify the ex situ versus in situ origin of the hot inner stellar halo by analyzing all TNG50 galaxies in the selection. We show the stellar ex situ fraction of the whole galaxy, f_ExSitu, versus the ex situ fraction of hot inner stellar halo, f_{ExSitu, halo}. The hot inner stellar halos in general have a higher ex situ fraction than the whole galaxy, although the two values become close in the galaxies with higher (global) ex situ fractions. However, in most galaxies, particularly those with lower masses (i.e., lower ex situ fractions), the ex situ fraction of stars in the hot inner stellar halo is smaller than 50%⁵.

All the findings presented above explain why a correlation naturally emerges between the mass of the most massive merger and the mass of the hot inner stellar halo of a galaxy, why such a correlation is stronger the higher the galaxy’s stellar mass, and why, at the same time, the relationship exhibits a lingering non-negligible galaxy-to-galaxy variation. Determining the origin of such a variation is the next and final objective of the paper.

5.4. Scatter of the correlations

5.4.1. Ex situ fraction of the hot inner stellar halo

In Fig. 7 we show that there is a strong correlation between the mass of the hot inner stellar halo, M_{*, halo(r < 2 Re)}, and the total ex situ stellar mass, M_*, Exsitu. However, there is still scatter, especially in galaxies below the average mass-size relation (g4, g7, g8). This scatter could be related to the ex situ fraction of galaxies: how many stars are accreted, how much the galactic structures are altered, and how much star formation is triggered by the merger events. As we have already seen in Fig. 5, the outliers are mostly galaxies with f_ExSitu < 0.1, and the majority of galaxies with different ex situ fractions are closely distributed along the correlations.

In order to understand the scatter, we separated the hot inner stellar halo into ex situ and in situ stars. In Fig. 13 we show the correlation between the total ex situ stellar mass, M_*, ExSitu, with M_{*, halo, ExSitu} (top) and M_{*, halo, InSitu} (bottom), respectively, for galaxies in g6 and g7. As a reminder, we note that g6 is the group of galaxies with large masses (10¹¹ < M_* < 10^11.6), average sizes, R_e, and with a strong correlation of M_{*, halo(r < 2 Re)} versus M_*, ExSitu, as shown in Fig. 7, while g7 is the group with small masses (10^10.3 < M_* < 10^10.6), small sizes (R_e below the average), and with the weakest correlation of M_{*, halo(r < 2 Re)} versus M_*, ExSitu.

Fig. 13. Dependence of the relation between M*, halo(r < 2 Re) and M*, ExSitu on the ex situ fraction of the hot inner stellar halo, fExSitu, halo. We split the hot inner stellar halo into two parts: M*, halo, ExSitu in the top panels and M*, halo, InSitu in the bottom. Here we show results for TNG50 in subsamples g6 and g7, color coded by fExSitu, halo. The plus symbols mark the galaxies with fExSitu, halo < 0.1. In g6, both M*, halo, ExSitu and M*, halo, InSitu are tightly correlated with the total ex situ stellar mass. In g7, M*, halo, ExSitu shows a similarly strong correlation with the total ex situ stellar mass; however, the trend with M*, halo, InSitu is weak with large scatter caused by the galaxies with low ex situ fractions of the hot inner stellar halo.

The data are color coded by the ex situ fraction of stars in the hot inner stellar halo, f_{ExSitu, halo}, and not by the f_ExSitu of the galaxy.

In g6, both M_{*, halo, ExSitu} and M_{*, halo, InSitu} are tightly correlated with the total ex situ stellar mass, although the correlation is slightly weaker for the latter. The combination of these two components leads to the strong correlation of M_{*, halo(r < 2 Re)} versus M_*, ExSitu. For these galaxies, mergers play the dominant role for the origin of both ex situ and in situ stars in the hot inner stellar halo. As we show in Sect. 5.2, in situ stars in the hot inner stellar halo should mainly be induced by mergers from a preexisting bulge and/or disk and formed during the star formation triggered by mergers.

In g7, M_{*, halo, ExSitu} shows a similarly strong correlation with the total ex situ stellar mass. However, the trend with M_{*, halo, InSitu} is weak. There are some galaxies with f_{ExSitu, halo} < 0.1, and, as we discuss in Sect. 3.2.2 and show in Fig. B.1, there is no clear transition from bulge to hot inner stellar halo at r ∼ 3.5 kpc for galaxies with low ex situ fractions (i.e., with f_ExSitu < 0.1, including those with f_{ExSitu, halo} < 0.1). For these galaxies, the origin of the in situ stars in the hot inner stellar halo is largely independent of mergers and instead due to the intrinsic scattering of the density distribution of the bulge. As shown in the bottom-right panel of Fig. 13, all galaxies with f_{ExSitu, halo} < 0.1 lie below the correlation and cause the large scatter in M_{*, halo, InSitu} versus M_*, ExSitu. The M_{*, halo, InSitu} is dominating the mass of the hot inner stellar halo, thus further causing the large scatter in the relation of M_{*, halo(r < 2 Re)} versus M_*, ExSitu.

5.4.2. Dominance of the most massive merger

In both the Milky Way and M 31, the buildup of the inner stellar halo is thought to be dominated by one massive merger event (Deason et al. 2015; Helmi et al. 2018; D’Souza & Bell 2018a). In more massive galaxies, the stellar halos are more likely to be built up through several large, equally massive mergers (Cooper et al. 2013). In our sample, which spans a wide mass range, the scatter of hot inner stellar halo mass, M_{*, halo(r < 2 Re)}, versus the mass of the most massive merger, M_*, Ex1, is further caused by the diversity of accretion histories and by whether and how much star formation is triggered during mergers. Here we show that the scatter is driven to a great extent by the number of massive merger events that have occurred over time.

In Fig. 14 we again show this correlation for galaxies in subsamples g6 and g7, but this time color coded by the contribution of the most massive merger to the total ex situ stellar mass: M_*, Ex1/M_*, ExSitu. Galaxies with accretion histories dominated by one massive merger tend to lie above the correlation: almost all the galaxies with M_*, Ex1/M_*, ExSitu > 0.8 lie above the median line. This is the case of TNG50 468590, shown as one case study (see Fig. 10). On the other hand, galaxies with multiple massive mergers tend to lie below the median of the correlation, such as TNG50 117253 (Fig. 11).

Fig. 14. Dependence of the relation between M*, halo(r < 2 Re) and M*, Ex1 on the dominance of the most massive merger event. Here we show results for TNG50 in subsamples g6 and g7, color coded by the ratio between the ex situ stellar mass brought in by the most massive merger and the total ex situ stellar mass: M*, Ex1/M*, ExSitu. The larger this ratio is, the larger the contribution of the most massive merger to the formation of the hot inner stellar halo.

6. Summary

In this paper we have analyzed the cosmological galaxy simulations TNG50, TNG100, TNG300, and EAGLE at z = 0. We find that the hot inner stellar halo of a galaxy, defined by the stars on dynamically hot orbits with 3.5 kpc < r < 2 R_e and λ_z < 0.5, is a strong predictor of the total ever-accreted stellar mass and of the stellar mass accreted from the most massive merger the galaxy has ever experienced.

In particular, we find a strong correlation between the hot inner stellar halo mass of a galaxy, M_{*, halo(r < 2 Re)}, and the total ex situ (i.e., accreted) stellar mass, M_ExSitu. The analog correlation of M_{*, halo(r < 2 Re)} versus the stellar mass of the most massive merger the galaxy has ever accreted, M_Ex1, is also strong.

These correlations hold almost unchanged in simulated galaxies throughout the inspected galaxy mass range from TNG50 and TNG100 (M_* in the range of 10^{10.3 − 11.6} M_⊙), TNG300 (at M_* > 10¹¹ M_⊙), and EAGLE (at M_* > 10^10.6 M_⊙). Therefore, our results are robust against different galaxy formation (i.e., feedback) models and across varying numerical resolution. The main result holds for those galaxy formation models and resolution choices that reproduce to a reasonable level of degree the most fundamental observed galaxy properties and statistics, as is the case for the TNG and EAGLE simulations. We find that TNG300 galaxies at stellar masses M_* ≲ 10¹¹ M_⊙ and EAGLE galaxies at M_* ≲ 10^10.6 M_⊙ deviate toward larger M_{*, halo(r < 2 Re)} at fixed satellite mass, and this is likely caused by an overestimation of M_{*, halo(r < 2 Re)} due to their larger bulge sizes.

Using the outcome of TNG50, we also find that the correlation between the inner-halo mass and the total ex situ stellar mass (or mass of the most massive merger) is strongest for massive and extended galaxies with M_* > 10^10.6 M_⊙ and R_e > 3 kpc, with scatter on the order of 0.1−0.2 dex; this scatter increases by an additional 0.1−0.3 dex in galaxies of lower mass and smaller stellar size.

The analysis of TNG50 galaxies clearly shows that the hot inner stellar halo is a product of the most massive merger(s) a galaxy has ever experienced. For galaxies with low ex situ stellar mass fractions, the radial probability density distributions of stars on dynamically hot orbits are highly concentrated within ∼3.5 kpc; on the other hand, galaxy mergers, and hence higher ex situ stellar mass fractions, produce an extended component that contributes to the hot inner stellar halo. We find that, according to TNG50, the stars in the hot inner stellar halo can have three main origins:

1. They can be accreted from satellite and merging galaxies.

2. They can be scattered out of the bulge and/or disk of a galaxy because of the interaction with satellite and merging galaxies.

3. They can be formed during the merger events and interactions.

The ex situ fraction within the hot inner stellar halo, origin pathway (1) of the inner-halo stars, is ∼50% in the most massive galaxies but can be smaller in lower-mass galaxies. Hence, our analysis shows that the formation of the hot inner stellar halo is not only due to the accretion of ex situ stars from merging and satellite galaxies but is also driven by the dynamical processes that merger events can trigger, such as enhanced star formation episodes and displacement and heating of in situ stars.

The mass of the hot inner stellar halo of galaxies defined in this paper is a quantity that we can robustly obtain from observations, as we will show in a companion paper. In fact, the maximum extension we advocate to define the hot inner stellar halo in this paper, 2 R_e, is set by limitations of observational data coverage rather than physical motivations. Using the correlations quantified here we will show in a companion paper how we can infer the stellar mass of the most massive merger(s) ever accreted by galaxies in the nearby Universe.

Acknowledgments

We thank the useful discussion with people from the Fornax 3D team: Enrica Iodice, Francesca Pinna, Ignacio Martin Navarro, Lodovico Coccato, Jesús Falcón-Barroso, Enrico Maria Corsini, Dimitri A. Gadotti, Katja Fahrion, Mariya Lyubenova, Richard McDermid, Adriano Poci, Mark Sarzi, and Tim de Zeeuw. L. Zhu 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. A.P. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subproject A01). G.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). D.N. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group (grant number NE 2441/1-1). 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: Comparison of time-averaged and phase-space-averaged stellar orbit distributions

In Fig. A.1 we show the stellar orbit distribution, p(r, λ_z), calculated in three different ways for an example galaxy, TNG50 468590: (1) from the instantaneous positions and velocities of particles at z = 0, (2) from phase-space-averaged values, and (3) from values that are time-averaged via the integration of the orbits in the frozen potential. The hot inner stellar halo fractions calculated from the three panels are 21.5%, 24.9%, and 24.5%, respectively, and the last two agree well with each other. This suggests that, for the purposes of the analysis here, it is sufficient to characterize galaxies in the (r, λ_z) plane via phase-space averaging.

Fig. A.1. Stellar orbit distribution, p(r, λz), calculated in three different ways for an example galaxy from TNG50 at z = 0. In the left, r and λz are calculated from the instantaneous positions and velocities of particles at z = 0. In the middle, r and λz are phase-space-averaged. In the right, r and λz are time-averaged by integrating the orbits in the frozen potential. Most of the orbits in the bulge and in the hot inner stellar halo are on box orbits, and particles on the box orbits span a wide range of λz but with time-averaged values of zero. By doing a phase-space averaging, we can narrow down the λz distribution of particles on box orbits but cannot shrink them to zero.

Appendix B: Distribution of r|_{max(CDF deviation)}

We can construct the deviation of CDF_hot from a control model, CDF_{hot, control}, for all TNG50 galaxies from subsamples g1-g9. In each group, CDF_{hot, control} is constructed by galaxies with f_ExSitu < 0.05 + min(f_ExSitu) in that group, in a similar way as shown in Fig. 3 for g6. As described in the main text of the paper, we define the radius where CDF_hot has the maximum deviation compared to the control model as r|_{max(CDF deviation)}, which reflects the transition from “bulge” to the merger-induced hot inner stellar halo for galaxies with relatively high ex situ fractions. In Fig. B.1 we show the distribution of r|_{max(CDF deviation)} for TNG50 galaxies.

Fig. B.1. Distribution of r|max(CDF deviation), defined as the radius where CDFhot − CDFhot, control reaches the maximum, for TNG50 galaxies in g1-g9. The black histogram represents galaxies with higher ex situ fractions (fExSitu > 0.1), and the red histogram galaxies with low ex situ fraction (fExSitu < 0.1).

For galaxies with f_ExSitu > 0.1, r|_{max(CDF deviation)} has a narrow distribution around 3.5 kpc: this is the case for galaxies in subsamples g2, g5, g6, g9, and, with larger variations, g3. The r|_{max(CDF deviation)} is smaller and distributed around 2.5 kpc for galaxies in g1 and g4 and around 2 kpc in g7. For galaxies with f_ExSitu < 0.1, r|_{max(CDF deviation)} has a wide distribution in all groups of galaxies, and it cannot reflect the real bulge size. In these galaxies, the variation in the density distribution of hot orbits is dominated by the intrinsic variation in the bulge and is not induced by mergers.

Appendix C: Correlation between the total stellar mass and accreted satellite mass

As a reference point to the main result of this paper (Fig. 5), we show the correlation of a galaxy’s stellar mass, M_*, with the stellar mass accreted from the most massive mergers, M_*, Ex1, for all TNG50 galaxies with 10^10.3 ≤ M_* ≤ 10^11.6 M_⊙ in the main panel of Fig. C.1. In comparison to the correlation between M_{*, halo(r < 2Re)} and M_*, Ex1 (Fig. 5, right panel), the correlation with total galaxy mass is weaker: the correlation coefficient reads ℛ = 0.75 and the 1σ scatter is about 0.1 dex larger in general. For completeness, in the upper panel of Fig. C.1 we also show the typical fraction of stellar mass contained within the hot inner stellar halos to the total stellar mass, according to our definition of hot inner stellar halo and to TNG50.

Fig. C.1. Correlation between a galaxy’s stellar mass, M*, and the stellar mass accreted from the most massive merger according to TNG50. Panel (a): Mass fraction of the hot inner stellar halo, M*, halo(r < 2Re), compared to the total stellar mass, M*, as a function of M*. The solid line is the average, and the dashed lines are ±1σ scatter. In this mass range, the hot inner stellar halo mass contributes only ∼5%−30% of the total stellar mass. Panel (b): Similar to the right panel of Fig. 5, but using the total stellar mass, M*, instead of the hot inner stellar halo mass along the x axis. Compared to Fig. 5, the correlation becomes weaker with ℛ = 0.75 and the 1σ scatter is about 0.1 dex larger in general. The ex situ fraction does not increase as gradually with increasing x as in Fig. 5, and the correlation of M* with fExSitu is also weaker than M*, halo(r < 2Re).

All Tables

All Figures

Fig. 1. Stellar mass and size of the simulated galaxies analyzed in this paper. Left and central panels: total stellar mass, M*, versus projected half-light size in the z band, Re, z band, of TNG50 galaxies at z = 0. The dashed black curve is the median size of TNG50 galaxies as a function of M*. The red and blue curves show the median sizes of early-type and late-type SDSS galaxies (Shen et al. 2003) for reference. TNG50 galaxies with 1010.3 < M* < 1011.6 M⊙ are color coded by their stellar ex situ fraction (see also the zoomed-in version in the middle panel). The ex situ fraction systematically increases with increasing stellar mass and size. However, there are still large variations among galaxies of similar mass and size. We further divide the TNG50 samples into nine subsamples as indicated by the dashed squares labeled g1–g9 from left to right and from top to bottom. Right panel: median stellar sizes as a function of M* for TNG50, TNG100, TNG300, EAGLE, and Illustris galaxies at z = 0. Here we use the 3D half-mass size, r1/2,  3D,  mass, for all galaxies. The vertical dashed lines mark the mass regions we focus on in this paper.

In the text

Fig. 2. Definition of hot inner stellar halo advocated for and adopted in this paper. We show the stellar orbit distribution, p(r, λz), stellar age distribution, t(r, λz), and metallicity distribution, Z(r, λz), in the phase-space of circularity, λz, versus radius, r, of a simulated galaxy, TNG50 468590. We define as a hot inner stellar halo the galactic component whose stars have orbits with rcut < r < rmax and λz < λz, cut, where rcut, rmax, and λz, cut are, respectively, 3.5 kpc, max(7 kpc, 2 Re), and 0.5.

In the text

Fig. 3. Properties of the stellar orbit distributions of simulated galaxies and physical justification of the definition of hot inner stellar halo. Top: probability density distribution in panel a and the average metallicity in panel b, as a function of radius, r, obtained by summing up all the hot orbits with λz < 0.5 from a selection of TNG50 galaxies at z = 0, namely those in group g6 of Fig. 1. Each galaxy’s curve is color coded by the galaxy’s stellar ex situ fraction, as indicated by the colorbar. The vertical line indicates the radius rcut = 3.5 kpc, where both the probability density and metallicity profiles flatten. Middle: distributions of the radii where the probability density distribution of hot orbits peaks, r|max(phot), for galaxies with M* > 1010.3 M⊙ in the TNG50, TNG100, TNG300, EAGLE, and Illustris simulations. Most galaxies in TNG50 have r|max(phot) < 1.5 kpc, consistent with real galaxies from the Fornax3D survey. Galaxies in TNG100, TNG300, and EAGLE have slightly larger bulge sizes but still with r|max(phot) < 3.5 kpc. Bottom: panel d: cumulative distribution function of the hot orbits with λz < 0.5, CDFhot as a function of radius, r, for a selection of TNG50 galaxies. Symbols and selection are the same as in the top panels. The green curve represents the control model constructed by the galaxies with the lowest stellar ex situ fraction. Bottom: panel e: deviation of CDFhot from the control model, CDFhot, control. The vertical line indicates the radius rcut = 3.5 kpc, which we adopt throughout as the inner boundary of the hot inner stellar halo.

In the text

Fig. 4. Connection between the hot inner stellar halo and the ex situ mass fraction according to TNG50. We show the mass fraction of the hot inner stellar halo compared to the total stellar mass within the maximum radius of rmax = 2 Re (fhalo = M*, halo(< 2 Re)/M*(< 2 Re)) versus the stellar ex situ fraction (fExSitu = M*, ExSitu/M*) for galaxies with M* > 1010.3 of TNG50 at z = 0. Here we exclude 37 objects identified with ongoing mergers or those with counter-rotating thin disks. Galaxies with size Re below, around, and above the average are indicated by dots, asterisks, and crosses, respectively: they are color coded based on stellar mass, with increasing stellar mass from black, to blue, to red. Galaxies with larger stellar masses and sizes exhibit systematically larger ex situ fractions as well as larger mass fractions in the hot inner stellar halo. The Pearson correlation coefficient ℛ(fhalo, fExSitu) = 0.86 of the whole sample is annotated. The black line is a linear fit to the data points, while the 1σ scatter of the data points against the linear fit is σ = 0.07.

In the text

Fig. 5. Correlations between the mass of the hot inner stellar halo, M*, halo(r < 2 Re), and the total accreted stellar masses (left) and the stellar mass accreted from the most massive merger (right), according to TNG50. Left panel: M*, ExSitu includes all the ex situ stellar mass ever accreted from all past mergers and from the stripping of ongoing mergers, satellites, and flybys. Right panel: M*, ExSitu is the stellar mass accreted from the most massive merger. Here all galaxies with M* > 1010.3 from TNG50 at z = 0 are included, color coded by their stellar ex situ fraction. The thick black curve is the running median, and the thin dashed curves represent the ±1σ scatter. Black plus symbols mark galaxies with ongoing merging or those with counter-rotating disks. The Pearson correlation coefficients, ℛ, of the correlations are annotated in each panel. The dashed magenta line refers to 1:1. Some satellite mergers will deposit their stellar mass across a wide radius, with only a fraction of it within 2 Re, and so can lie above that line. On the other hand, the x axis shows a combination of in situ and ex situ stars such that the total mass can greatly exceed the total mass of the satellites in some galaxies with low ex situ fractions, causing them to lie below the line.

In the text

Fig. 6. Information content of the hot inner stellar halo across galaxy formation simulations. We show the correlations of M*, halo(r < 2 Re) with the total ex situ stellar mass, M*, ExSitu (left), and with the stellar mass accreted from the most massive merger, M*, Ex1 (right), for different galaxy simulations. The thick black, blue, yellow, and green lines are the median curves from TNG50, TNG100, TNG300, and EAGLE galaxies at z = 0, respectively, in the M* = 1010.3 − 11.6 M⊙ range; the thin curves are the ±1σ scatters. The Pearson correlation coefficients, R, from the four sets of simulations are labeled in the corresponding colors.

In the text

Fig. 7. Correlation between M*, halo(r < 2 Re) and the total ex situ stellar mass, M*, ExSitu, as a function of galaxy stellar mass and stellar size. Here we show the relationships for TNG50 galaxies divided into nine groups by mass and size, as shown in Fig. 1, one grouping per panel, with galaxy stellar mass increasing from left to right and galaxy stellar size increasing from bottom to top. In the middle row, galaxies of median size for their mass are depicted. The black dots denote TNG50 galaxies without ongoing mergers, whereas the gray plus symbols indicate ongoing mergers. The thick black line is a linear fit to the black dots, and the thin black lines represent the ±1σ scatter. The median and ±1σ scatter of the whole sample as shown in Fig. 5 are reproduced as gray curves for reference. The medium and ±1σ uncertainty of the Pearson correlation coefficient, ℛ(M*, halo(r < 2 Re), M*, ExSitu), slope of the linear fit, p, and the typical 1σ scatter of the data points against the linear fit are labeled for each group. The four galaxies that are discussed as case studies in Sect. 5 are marked with colored symbols (in panels g6 and g9).

In the text

	Fig. 8. Correlation between the hot inner stellar halo mass, M*, halo(r < 2 Re), and the stellar mass accreted from the most massive merger, M*, Ex1, as a function of galaxy mass and size. Annotations are as in Fig. 7, and results are from TNG50.
In the text

Fig. 9. Dependence of the correlation between the masses of the most massive merger and of the hot inner stellar halo of a galaxy on the operational definition of the latter. We show the Pearson correlation coefficient between M*, Ex1 and M*, halo(r < 2 Re) as a function of λz, cut, the maximum level of orbital “coldness” adopted to define inner-halo stars. Namely, in the plot we adopt different definitions of hot inner stellar halo by imposing λz < λz, cut. Dashed curves in different colors indicate galaxies in subsamples g1–g9, and the solid black curve indicates galaxies of the whole sample.

In the text

Fig. 10. Formation of the hot inner stellar halo in TNG50 468590 (a-1 and b-1; TNG50 468590 is the same galaxy as illustrated in Fig. 2) and TNG50 375073 (a-2 and b-2). Panel a-1: stellar mass assembly history of TNG50 468590. The black dots indicate the stellar mass evolution of the main progenitor as a function of redshift, z. The dashed colored lines denote the mass assembly of the secondary galaxies that merged into the main progenitors directly and at one point reached a maximum stellar mass larger than 109 M⊙; the corresponding colored star symbols indicate the times when the maximum mass was reached. TNG50 468590 experienced a massive merger at z ∼ 1 with a stellar mass ratio of ∼1:1. The solid and dashed red vertical lines indicate the beginning and ending of this merger. Top row of panel b-1: morphology of the main progenitor before the major merger in the leftmost panel and the morphology at redshift z = 0 in the rightmost panel; in between we have the contribution to the z = 0 morphology from different origins with (from left to right) the stars from the main progenitor, accreted from the most massive merger, formed during the major merger, formed after the major merger ends, and accreted from subsequent minor mergers. We use co-moving distances in units of kpc (labeled as ckpc) for galaxies at z > 0. The surface density is in units of stellar mass per unit area kpc2. The stellar mass is normalized such that the total stellar mass of the galaxy at z = 0 equals unity within the figure coverage. The bottom row shows the orbital distribution of the corresponding stars in the phase-space of circularity, λz, versus radius, r, as derived from the simulation particle information. The probability density is in units of stellar mass per unit area in the phase-space. The dashed lines indicate our orbital-based division into four components as adopted for NGC 1380: disk, warm component, bulge, and hot inner stellar halo. It has a hot inner stellar halo mass of 3.4 × 1010 M⊙, and 90% of it was produced by the most massive merger event. The contribution from subsequent minor mergers is negligible. Panels a-2 and b-2: are similar but for TNG50 375073. TNG50 375073 had two massive mergers, one at z ∼ 2 and one at z ∼ 1. Panel b-2: we trace the structure formation during the most massive merger, which happened at z ∼ 1. It has a hot inner stellar halo mass of 3.0 × 1010 M⊙, and 84% of it was produced by the most massive merger event.

In the text

Fig. 11. Formation of the hot inner stellar halo in the simulated galaxies TNG50 2 (a-3 and b-3) and TNG50 117253 (a-4 and b-4). The symbols are the same as in Fig. 10. Panel a-3: TNG50 2 has accreted three massive satellites, which are denoted by Ex1, Ex2, and Ex3, with decreasing stellar mass. Here, Ex3 and Ex2 were accreted simultaneously at z ∼ 1 and marked as the first merger, and Ex3 was accreted at z ∼ 0.1 and marked as the second merger. Panel b-3: we trace the structure formation during the two merger events. During the first merger, stars accreted from both Ex2 and Ex3, as well as stars formed during the merger, contribute to the hot inner stellar halo. After the first merger and before the second merger, a large number of stars were formed in a cold disk. The disk was then heated, becoming warm, during the second merger, and stars accreted from Ex1 contribute significantly to the hot inner stellar halo. Stars are rarely formed during and after the second merger as there is little gas left. TNG50 2 has a hot inner stellar halo mass of 4.3 × 1010 M⊙, of which 45% is a product of the most massive merger. Panel a-4: TNG50 117253 has accreted three massive satellites, which are still denoted by Ex1, Ex2, and Ex3, with decreasing stellar mass. We focus on the last two accretion events, Ex1 and Ex2, at 1 < z < 1.8. Panel b-4: structure formation during the accretion of Ex1 and Ex2. TNG50 117253 has a hot inner stellar halo mass of 1.5 × 1010 M⊙: the fraction produced by the most massive merger has an upper limit of 60%.

In the text

Fig. 12. Ex situ and in situ origin of the hot inner stellar halo of TNG50 galaxies. Top: probability density distribution of stars on hot orbits, for in situ (left) and ex situ (right) stars separately. Here we divide the subsample of TNG50 galaxies in g6 into five bins according to their ex situ fraction, fExSitu. Each curve represents the average phot (in situ or ex situ) of galaxies in each bin, color coded by its range of stellar ex situ fraction from 0−0.2 to 0.5−0.7, as indicated by the colored numbers. The curves are normalized such that the total probability of all stars within 2 Re equals unity. Bottom: stellar ex situ fraction of the galaxy, fExSitu, versus the ex situ fraction of stars in the hot inner stellar halo, fExSitu, halo (3.5 kpc < r < 2 Re and λz < 0.5). Each dot represents one galaxy in TNG50, color coded by galaxies stellar mass, M*, and crosses denote galaxies in subsample g6 (i.e., with a galaxy stellar mass in the range 1011 − 1011.6 and an average stellar size).

In the text

Fig. 13. Dependence of the relation between M*, halo(r < 2 Re) and M*, ExSitu on the ex situ fraction of the hot inner stellar halo, fExSitu, halo. We split the hot inner stellar halo into two parts: M*, halo, ExSitu in the top panels and M*, halo, InSitu in the bottom. Here we show results for TNG50 in subsamples g6 and g7, color coded by fExSitu, halo. The plus symbols mark the galaxies with fExSitu, halo < 0.1. In g6, both M*, halo, ExSitu and M*, halo, InSitu are tightly correlated with the total ex situ stellar mass. In g7, M*, halo, ExSitu shows a similarly strong correlation with the total ex situ stellar mass; however, the trend with M*, halo, InSitu is weak with large scatter caused by the galaxies with low ex situ fractions of the hot inner stellar halo.

In the text

Fig. 14. Dependence of the relation between M*, halo(r < 2 Re) and M*, Ex1 on the dominance of the most massive merger event. Here we show results for TNG50 in subsamples g6 and g7, color coded by the ratio between the ex situ stellar mass brought in by the most massive merger and the total ex situ stellar mass: M*, Ex1/M*, ExSitu. The larger this ratio is, the larger the contribution of the most massive merger to the formation of the hot inner stellar halo.

In the text

Fig. A.1. Stellar orbit distribution, p(r, λz), calculated in three different ways for an example galaxy from TNG50 at z = 0. In the left, r and λz are calculated from the instantaneous positions and velocities of particles at z = 0. In the middle, r and λz are phase-space-averaged. In the right, r and λz are time-averaged by integrating the orbits in the frozen potential. Most of the orbits in the bulge and in the hot inner stellar halo are on box orbits, and particles on the box orbits span a wide range of λz but with time-averaged values of zero. By doing a phase-space averaging, we can narrow down the λz distribution of particles on box orbits but cannot shrink them to zero.

In the text

	Fig. B.1. Distribution of r|max(CDF deviation), defined as the radius where CDFhot − CDFhot, control reaches the maximum, for TNG50 galaxies in g1-g9. The black histogram represents galaxies with higher ex situ fractions (fExSitu > 0.1), and the red histogram galaxies with low ex situ fraction (fExSitu < 0.1).
In the text

Fig. C.1. Correlation between a galaxy’s stellar mass, M*, and the stellar mass accreted from the most massive merger according to TNG50. Panel (a): Mass fraction of the hot inner stellar halo, M*, halo(r < 2Re), compared to the total stellar mass, M*, as a function of M*. The solid line is the average, and the dashed lines are ±1σ scatter. In this mass range, the hot inner stellar halo mass contributes only ∼5%−30% of the total stellar mass. Panel (b): Similar to the right panel of Fig. 5, but using the total stellar mass, M*, instead of the hot inner stellar halo mass along the x axis. Compared to Fig. 5, the correlation becomes weaker with ℛ = 0.75 and the 1σ scatter is about 0.1 dex larger in general. The ex situ fraction does not increase as gradually with increasing x as in Fig. 5, and the correlation of M* with fExSitu is also weaker than M*, halo(r < 2Re).

In the text