Discovery of an ancient massive merger event in the Fornax cluster galaxy NGC 1380

in the galaxy 1380 in the a population-orbital superposition metallicity VLT/MUSE and metallicity distributions of an inner with metallicity Z/Z

Driven by gravity, galaxies are expected to continuously grow through the merging of smaller systems. To derive their past merger history is challenging, as the accreted stars disperse quickly; yet, it is a needed step to test the theory of hierarchical evolution. The merger histories of the most massive Local Group spirals, the Milky Way and M31, have been recently uncovered by using the motion and chemistry of their individual stars. On the other hand, the details of the merger history of galaxies at larger distance have so far remained hidden. Here we report the discovery of an ancient, massive merger event in the lenticular galaxy NGC 1380 in the Fornax cluster. By applying a recently developed populationorbital superposition model 1 to NGC 1380's surface brightness as well as stellar kinematic, age, and metallicity maps from VLT/MUSE IFU data 2 , we obtain the stellar orbits, age and metallicity distributions of this galaxy. The highly radial orbits which make up an inner stellar halo are ∼ 13 Gyr old with metallicity Z/Z ∼ 1.2 and comprise a stellar mass of M * ,halo(r<2Re) ∼ 3.4 × 10 10 M . By comparing to analogues from the cosmological galaxy simulation TNG50 3,4 , we find that the formation of the inner stellar halo of NGC 1380 requires a merger with a massive satellite galaxy with stellar mass of ∼ 3 × 10 10 M that occurred roughly ∼ 10 Gyr ago. Moreover, we infer the total accreted stellar mass of NGC 1380 to be ∼ 6 × 10 10 M . The massive merger in NGC 1380 is the first major merger event found in a normal phase-mixed galaxy beyond the Local Volume, and it is the oldest and most massive one identified in nearby galaxies so far. Our chemo-dynamical method, when applied to extended deep IFU data and in combination with cosmological galaxy simulations, can quantitatively unravel the merger history of a large number of nearby galaxies.
NGC 1380 data and model. NGC 1380 (FCC 167) is a lenticular (S0) galaxy in the Fornax cluster at a distance of 21.2 Mpc 5 , so that on the sky 1 102 pc. Deep photometry from the Fornax Deep Survey 6 with the VLT Survey Telescope reveals a total magnitude in the r band of m r = 9.27 and half-light-radius of R e = 56.2 . The total stellar mass inferred from the bestfit orbit-based dynamical model (see Methods) is M * = 1.8 × 10 11 M . NGC 1380 is observed with the VLT Multi Unit Spectroscopic Explorer (MUSE) instrument as part of the Fornax3D survey 2 . From these deep integral-field unit (IFU) data, we extract high-quality stellar kinematic maps as well as light-weighted stellar age and metallicity maps. We then create population-orbit superposition models and fit those to the observed surface brightness of NGC 1380 as well as the extracted stellar kinematic, age and metallicity maps (see Methods). Our best-fitting models match very well the stellar population and kinematics data as shown in Figure 1.
Our method yields the weights of the different stellar orbits that contribute to the best-fitting model as well as the average stellar age and metallicity of the stars represented by each orbit. We next characterize each stellar orbit with two main properties 7 : the time-averaged radius, r, which represents the size of the orbit, and the circularity, λ z = J z /J max (E), which represents the angular momentum of the orbit around the minor z-axis normalized by the maximum of a circular orbit with the same binding energy E. While |λ z | ∼ 1 represents dynamically cold orbits dominated by regular rotation, |λ z | ∼ 0 represents dynamically hot orbits which are dominated by radial random motions or long axis tubes with regular rotation around the long axis. Negative values λ z < 0 refer to orbits that counter-rotate with respect to the net (prograde) rotation. Figure 2 presents the resulting stellar orbit distribution of NGC 1380 as a probability density distribution p(r, λ z ) in the phase-space r versus λ z , as well as the age t(r, λ z ) and metallicity Z(r, λ z ) distributions of the stellar orbits in the same phase-space.
Orbit-based decomposition of NGC 1380. Considering the sub-structures in the stellar orbit probability density, age and metallicity distributions in the phase-space of r versus λ z , we separate NGC 1380 into four components: disk (λ z > 0.8), bulge (λ z < 0.8, r < r cut ), warm (0.5 < λ z < 0.8, r > r cut ), and inner stellar halo (λ z < 0.5, r > r cut ). Here, r cut = 3.6 kpc 0.6R e for NGC 1380 is based on a transition in the density distribution of the dynamically hot orbits with an associated transition in metallicity (see Methods). The disk, bulge, warm, and inner stellar halo components contribute respectively 19, 38, 17 and 26 per cent to the total r-band luminosity of the galaxy within 2 R e . As we tested with mock data created from simulations, the luminosity fractions of such four components are well recovered by our model (see Methods). The best-fit global stellar mass-to-light ratio M/L r = 2.7 M /L ,r implies that the inner stellar halo has a stellar mass M * ,halo(r<2Re) = 3.4 × 10 10 M (see Methods).
The morphology, age and metallicity profiles of each component can be reconstructed from our recovered distribution of orbits. In the bottom panels of Figure 2, we show profiles along the major axis of NGC 1380 of the surface brightness, stellar age, and metallicity for each of the four components. The disk is radially extended and has a nearly exponential profile; it is relatively young and metal-rich. The bulge is concentrated, old and metal-rich. The inner stellar halo is radially extended, also approximately described by an exponential profile; it is old and metal-poor comparing to the other components of this galaxy. The warm component is radially less extended but as old and similarly metal-poor as the inner stellar halo. The four components are decomposed in such a way that they have different structural, kinematical, age, and metallicity distributions and thus likely have different physical origins.
Inner stellar halos induced by massive mergers. To connect the origin of the inner stellar halo with the assembly history of NGC 1380, we turn to galaxy simulations. We select galaxies with stellar mass M * and half-light-radius R e similar to NGC 1380 from the cosmological galaxy simulation TNG50 3, 4 , the highest resolution realization of the IllustrisTNG project 1 . This results in 50 simulated galaxies without on-going mergers at z = 0: they typically consist of four phase-space components similar to NGC 1380. We calculate the stellar mass fraction of the four components within 2R e of each simulated galaxy, and choose five of them that have similar (or slightly less) inner stellar halo masses as NGC 1380 for an in-depth study (see Figure 5 in Methods).
All the five NGC 1380-like simulated galaxies have experienced at least one massive merger with M * ,acc1 > 10 10 M during their life time, where M * ,acc1 denotes the total stellar mass accreted from the companion during the merger. Their merger event may have lasted for as long as ∼ 2 Gyr.
We show the assembly history and the formation of structures during and after the merger for one of the simulated galaxies of reference, denoted TNG50 468590, in Figure 3 and for the other four study cases in the Methods. TNG50 468590 experienced a massive merger with stellar mass ratio of ∼ 1 : 1 at redshift z ∼ 1. Before the merger, the main progenitor was dominated by a cold disk and a compact bulge with no substantial inner stellar halo. During the major merger, both the main progenitor and the satellite were destroyed, with some stars settling into the bulge while others forming part of the inner stellar halo. Both the main progenitor and the satellite were gas-rich. As such, 4 × 10 10 M of stars (comparable to the mass of the main progenitor) formed during the merger, most of which contributed to the bulge, while ∼ 15 per cent of them were distributed into the inner stellar halo. After the merger, new stars formed on dynamically cold disk-like orbits. The latter disk structure persists 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 inner stellar halo. The inner stellar halo mass of TNG50 468590 at z = 0 is 3.3 × 10 10 M : importantly, only ∼ 10 per cent of the inner halo stars was already there in the main progenitor before the massive merger started; ∼ 30 per cent was in the bulge or disk of the main progenitor and redistributed to the inner stellar halo during the merger; ∼ 40 per cent was brought in by the satellite galaxy, and 20 per cent formed during the merger and immediately distributed into the inner stellar halo regions. In total, about 90 per cent of the inner stellar halo of this simulated galaxy was the product of the massive merger event.
TNG50 468590 is somewhat a special case, in that its history is characterized by only one dominating massive merger throughout cosmic epochs. Furthermore, the merger has a stellar mass ratio close to unity, so that both the main progenitor and the satellite contribute significantly to the bulge as well as to the inner stellar halo. The other four simulated galaxies with similar inner stellar halo mass to TNG50 468590 but with a variety of accretion histories are described in the Methods. In all cases, the stars accreted from the most massive satellite(s) favour 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 8 .
With the 50 simulated galaxies from TNG50, we find a strong correlation between the mass of their inner stellar halo and that of their accreted satellites ( Figure 4). We infer a linear correlation between the inner stellar halo and the most-massive, ever-accreted satellite of M * ,sat1 /10 10 M = 0.67 × M * ,halo(r<2Re) /10 10 M + 0.09, with a Pearson correlation coefficient of R(M * ,halo(r<2Re) , M * ,sat1 ) = 0.74. Here, the satellite mass M * ,sat1 is defined as the maximum stellar mass it has ever reached before the merger ended 9 . As a merger could last for a few Gyr and as a merging galaxy could be continuously loosing and forming new stars during the event, the total stellar mass accreted from such a satellite during the merger, M * ,acc1 is a priori different from M * ,sat1 , and in fact it is typically larger. Nevertheless, a similar correlation holds, with M * ,acc1 /10 10 M = 0.95 × M * ,halo(r<2Re) /10 10 M + 0.23 with R(M * ,halo(r<2Re) , M * ,acc1 ) = 0.70. Accounting for the stellar mass of all accreted satellites makes the correlation stronger, which gives M * ,all,acc /10 10 M = 1.71×M * ,halo(r<2Re) /10 10 M +0.27 with R(M * ,halo(r<2Re) , M * ,all,acc ) = 0.83. The outcome of the TNG50 simulation reveals that the mass of the inner stellar halo is a good predictor for the mass of the most massive satellites or the total stellar mass ever accreted by a galaxy observed at z = 0. These correlations hold amid the variations in merger times and merging orbits of the satellites, as we have only placed constraints on stellar mass M * and size R e of the adopted simulated galaxy sample.
The mass and timing of NGC 1380's past major merger. The inner stellar halo mass within 2 R e can thus be used to infer the mass of the most massive satellite as well as of the total accreted or ex-situ stellar mass of an observed galaxy, in combination with the fact that the inner stellar halo mass can be well recovered by our population-orbit based model for real galaxies even long after the stars from the merger events are phase mixed (see Methods).
From the TNG50 scaling relations, we can infer that, with an inner stellar halo mass of M * ,halo(r<2Re) = 3.4 × 10 10 M , the most massive satellite galaxy that NGC 1380 ever accreted had a maximum stellar mass of M * ,sat1 = 2.4 ± 0.7 × 10 10 M , with the total stellar mass accreted from this satellite being as large as M * ,acc1 = 3.5 ± 1.1 × 10 10 M . We check that these results are in excellent agreement with independent estimates from analytic merger histories incorporating chemical constraints of the system, which imply M * ,acc1 = (3 ± 1) × 10 10 M (Leaman et al., in prep). From the TNG50 scaling relations, in a similar fashion we can infer that NGC 1380's total accreted mass is M * ,all,acc = (6.1 ± 1.3) × 10 10 M . Here, we take the 1σ scatter of the corresponding fitted line shown in Figure 4 as error of the inferred satellite(s) mass here. The satellite(s) mass inferred from these relations do not change with slightly different r cut used for inner stellar halo separation (see Methods). These ex-situ stars can be distributed throughout the galaxy's stellar body, i.e., far beyond the observed field-of-view. We hence estimate that the accreted stellar mass fraction of NGC 1380 across the whole galaxy is 6.1 × 10 10 M /1.8 × 10 11 M , i.e. about 35 per cent of the total stellar mass of the galaxy.
We can further infer the time of the last massive merger from the two-dimensional distribution of stellar age and orbital circularity. A dynamically cold disk is usually, at least partially, destroyed during a massive merger, resulting in stars of similar age as the surviving disk stars but on dynamically warmer orbits. Only after the merger and in case of sufficient gas presence, newly formed stars can re-grow a stellar disk: this structure can persist for a long time if not disrupted by subsequent massive merger events, while minor mergers only heat up the disk without markedly affecting the circularity distribution 10 .
We find a transition in the circularity λ z versus stellar age t plane, p(λ z , t) for the distribution of stars in our model of NGC 1380, as shown in panel (f) of Figure 4: for ages younger than t 10 Gyr, suddenly stars are only on cold orbits with λ z > 0.8, i.e., the fraction of disk-like orbits, f disk , approaches unity. Adopting 1 − f disk = 0.5 for the transition, we infer that the massive merger responsible for the build up of NCG 1380's inner stellar halo ended about 10 Gyr ago (see Methods). This approach is supported by what we find with the simulated galaxy analogues. In panel (d) of Figure 4, we show the probability distribution of the stellar orbits of the simulated galaxy TNG50 468590 in p(λ z , t) plane. There is a sharp transition of stars from nondisk dominated orbits to disk dominated orbits: knowing every aspect of the merger history of the simulated galaxies and via connection to the stellar orbit distribution, we can confidently associate this transition to the time when the last massive merger ended. This is also true for all the five case-study galaxies (see Methods). In panel (e), we show how the same distributions are returned by the best-fit model to mock data of TNG50 468590. The transition is recovered, although it is not as sharp as in the simulation data due to uncertainties of the stellar age of orbits in our model. NGC 1380 is located in the so-called north-south clump of the Fornax cluster 11,12 , which is comprised of a group of galaxies that are thought to have accreted into the potential well of the massive cluster core more than 8 Gyr ago 13 . Two other early-type galaxies in the same clump, NGC 1380A (FCC 177) and NGC 1381 (FCC 170), have also been suggested to have accreted a substantial fraction of stars ∼ 9 Gyr ago 14,15,16 . As galaxy-galaxy mergers in cluster environments are unlikely because of the high relative velocity of their members 17 , it is expected that these galaxies underwent their massive mergers before they fell into the Fornax cluster. The merger of NGC 1380 with its most massive satellite occurring ∼ 10 Gyr ago is consistent with this scenario, and supports the idea of pre-processing through mergers in groups before falling into the cluster.
Comparison with the Milky Way and the Andromeda galaxy. In the Milky Way (MW), an inner stellar halo on highly radial orbits, the Gaia-Enceladus, was identified, which is thought to have originated from a single satellite M * ,acc1 ∼ 3 − 6 × 10 8 M , around 10 Gyr ago 18, 19 . In the Andromeda galaxy (M31), an inner stellar halo with a giant stream has been identified 20 , which could be also induced by a single massive merger, called M32p, with M * ,acc1 ∼ 2 × 10 10 M occurred ∼ 2 Gyr ago 21, 22 . This merger event also explains M31's dynamically warm disk with a ratio of the ordered to random motion of V /σ ≤ 3 23 (corresponding to λ z ≤ 0.8). In both the MW and M31, the quantification of the satellite mass relies on the mass-metallicity relation, while the constraint on merger time relies on star formation histories obtained from resolved stellar population of individual stars. Similar studies are thus limited to the Local Group until robust star formation histories from integrated spectra are possible for the oldest age stars, or until 30-m telescopes come and allow to get resolved stellar populations for more distant galaxies.
Here, we have circumvented the latter limitations and inferred a specific and important aspect of the merging history of NGC 1380, which is about thirty times further away than M31, via a new method. The inner stellar halo we identified proves to be a strong predictor of the most massive merger the system experienced, captured by a simple linear correlation between M * ,acc1 versus M * ,halo(r<2Re) , where the latter is well determined from our model. Moreover, the merger time can be well constrained by the age distribution of stars on cold, disk-like orbits in comparison to those on dynamically warmer orbits, in a way similar to MW studies.
Among the most massive merger events in the only three systems that have been uncovered in the nearby Universe, NGC 1380 has the oldest and most massive one (see Methods). Unlike the MW and M31, NGC 1380 is an early-type galaxy without on-going star formation in the disk. NGC 1380 has a total stellar mass about 1.5 times of M31; its most massive merger mass of ∼ 3 × 10 10 M is also about 1.5 times of M31's merger mass. However, the merger in NGC 1380 happened much earlier than in M31. Given the subsequent evolution in the cluster environment, the high encounter velocities and long relaxation times should have resulted in a smoother structural distribution, and possibly also a steeper density slope for NGC 1380's halo compared to that of M31 24 .
Similar IFU data as we used for NGC 1380 have been obtained for hundreds of nearby galaxies through observations with instruments like VLT/MUSE. Our discovery thus opens a new way to quantify the merger history of a large number of nearby galaxies and so place important constraints on the assembly of galaxies as a function of environment and in a cosmological context.   The first column shows the data from the VLT/MUSE observations, the second column the best-fit population-orbit model, and the third column the residuals of data minus model divided by the observational errors. From top to bottom: surface brightness, mean velocity, velocity dispersion, Gauss-Hermite coefficients h 3 and h 4 , as well as light-weighted age and metallicity maps of the stars in NGC 1380. Note that kinematic maps were derived from spectra binned with a different scheme than the age and metallicity maps, so their bin numbers are not the same. The probability density distribution of the stellar orbits p(r, λ z ) in the phase-space of radius r versus circularity λ z in unit of stellar mass per unit area in the phase space. We show the model within 140 arcsecs = 12 kpc for this galaxy. The stellar mass are normalised such that they sum to unity within the data coverage. Panels (b) and (c): The stellar age t(r, λ z ) and metallicity Z(t, λ z ) distribution of the orbits in the same phase-space. All distributions p(r, λ z ), t(r, λ z ), and Z(t, λ z ) are averages of best-fitting models that fall within the 1σ confidence level of the model hyper-parameters (see Methods The merger starts  Figure 3: The formation of the inner stellar halo in a NGC 1380 analogue from the TNG50 simulation: TNG50 468590. Panel (a): the stellar mass assembly history of TNG50 468590, according to the SUBLINK merger-tree code 9 . The black dots indicate the stellar mass evolution of the main progenitor as a function of redshift, z. The colored dashed lines denote the mass assembly of those secondary galaxies that merged into the main progenitors directly and have once reached a maximum stellar mass larger than 10 9 M ; the corresponding colored stars indicate the times when such maximum mass was reached. TNG50 468590 experienced a massive merger at z ∼ 1 with stellar mass ratio of ∼ 1 : 1. The solid and dashed red vertical lines indicate the beginning and ending of this merger. Panel (b): in the top row, the leftmost panel shows the morphology of the main progenitor before the major merger and the rightmost one the morphology at redshift z = 0. We use co-moving distances in units of kpc (labeled as ckpc) for galaxies at z > 0. In between are 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 satellite, formed during the major merger, formed after the major merger ends, and accreted from subsequent minor mergers. The surface density is in units of stellar mass per unit area kpc 2 . The stellar mass is normalised such that the total stellar mass of the galaxy at z = 0 sum to unity within the figure coverage. The bottom row shows the orbital distribution of the corresponding stars in the phasespace of radius r versus circularity λ z as derived from the simulation particle information. The probability density is in unit 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 inner stellar halo. It has an inner stellar halo mass of 3.3 × 10 10 M , and 90 per cent of it was produced by the most massive merger. The contribution from subsequent minor mergers is negligible.

Methods
Stellar kinematics, age, and metallicity maps of NGC 1380. NGC 1380 was observed, as part of the Fornax Deep Survey, with the VLT Survey Telescope down to 27 mag arcsec −2 in r-band and out to a radius of 10 R e . Photometric analysis 6 yields a total magnitude of m r = 9.27 and a halflight-radius of R e = 56.2 = 5.8 kpc. NGC 1380 was observed, as part of the Fornax3D survey 2 , by the VLT/MUSE instrument with three pointings out to a galactocentric radius of 2.7 R e 150 along the major axis and up to 0.7 R e 40 along the minor axis. The spectra in the latter IFU data were spatially binned using a Voronoi tessellation 26 to reach a minimal signal-to-noise S/N = 100, resulting in a total of 4300 bins. Following earlier analysis of NGC 1380 2 , the stellar kinematics were extracted by applying pPXF full-spectral fitting 27 to the wavelength range 4750-5500Å. This yields high-quality maps of the stellar mean velocity V , velocity dispersion σ, and higher order velocity moments parameterised through the Gauss-Hermite coefficients h 3 and h 4 .
To also obtain high-quality maps of stellar age and metallicity, the spectra were re-binned to reach S/N = 200, resulting in a total of 974 bins. Each binned spectrum was fitted with regularized pPXF 28 using the MILES 29 library with a range of stellar ages, metallicities and alpha-abundances in [Mg/Fe] for a constant MW-like Kroupa 30 stellar initial mass function. Following earlier work 14 , for each bin we obtained a light-weighted stellar age t and metallicity Z/Z .
Population-orbit superposition model of NGC 1380. For this paper, we have performed a population-orbit superposition model to fit the surface brightness, stellar kinematic, age, and metallicity maps of NGC 1380. Our method was introduced and extensively verified 1 , we only briefly describe the model construction for NGC 1380 here.
The gravitational potential is assumed to be generated by the stellar mass plus dark matter halo. Also a central black hole of 10 8 M is added although it remains unresolved by the kinematic data and does not affect our results. To obtain the stellar mass distribution, we first fitted the rband image by a Multiple Gaussian Expansion (MGE) model 31 . Then by adopting a set of viewing angles (ϑ, ψ, φ), we de-projected toward a triaxial MGE model which represents the intrinsic stellar luminosity density. Finally, by multiplying with a constant stellar mass-to-light ratio M/L r , we arrived at the intrinsic stellar mass distribution. The effects of the constant stellar mass-to-light ratio to our results are discussed below.
The three viewing angles relate directly to three axis ratios (p, q, u) describing the intrinsic shape 32 . We leave the intrinsic intermediate-to-major and minor-to-major axis ratios (p, q) as free parameters, but fix u = 0.9999 so that the intrinsic major axis projects onto the projected major axis. While the latter is formally only valid for oblate axisymmetric systems, it does not affect the results for NGC 1380 which is at most mildly triaxial. The dark halo is throughout assumed to be spherical with a NFW 33 radial profile with two free parameters the halo (virial) mass M 200 defined as the total mass enclosed within the virial radius r 200 within which the average density is 200 times the critical density, and concentration c defined as the ratio between dark matter virial radius and the dark matter scale radius. In total we thus have five free so-called hyper-parameters: M/L r , p, q, M 200 and c.
For each set of hyper-parameters, we computed tens of thousands of orbits in the corresponding gravitational potential. The weights of the orbits in the resulting orbit library were then determined from fitting simultaneously the orbit-superposition to the intrinsic and projected luminosity density and stellar kinematic maps. Once we have explored the hyper-parameter space, we selected the best-fitting models with ∆χ 2 ≡ χ 2 −min(χ 2 ) < √ 2 × n GH × N obs , where n GH = 4 is the number of stellar kinematic moments and N obs = 4300 is the number of bins of each kinematic map. There are 55 models in total within this 1σ confidence level. Our dynamical model constrain the total stellar mass of 1.33 × 10 11 M within 2 R e , with M/L r = 2.7 ± 0.3 M /L . The mass out of the kinematic data coverage is not well constrained, we just take two times of the mass within R e as NGC 1380's total stellar mass, which gives M * = 1.82 × 10 11 M .
Next, for each of the 55 models, we applied Voronoi binning in the phase space of r versus λ z to group the orbits into ∼ 150 bundles with comparable weight. Then we tagged each bundle with a single stellar age t and metallicity Z and compared the superposition of the orbit bundles with the light-weighted age and metallicity maps extracted from the VLT/MUSE observations. We used Bayesian statistical analysis (Python package pymc3) to obtain age and metallicity of the orbital bundles. What results from the Bayesian analysis is a distribution of age/metallicity for each orbit bundle: their average values are taken as the best-fit solutions. The combination of the best-fit orbital weights with best-fit age and metallicity of orbital bundles yields population-orbital models of NGC 1380 that fit in detail the observed stellar kinematic, age, and metallicity maps.
Extensive testing on simulated galaxies 1 shows that our population-orbital models are able to robustly recover the internal stellar orbit distribution in r versus λ z , the stellar population distribution in t versus Z, and the correlation between t and λ z . When galaxies are viewed closer to face-on, the recovery is more difficult, but can be aided by a theoretically motivated age-metallicity relation. The latter is not needed in case of NGC 1380 which is close to edge-on at an inclination angle ϑ = 76 • . TNG50 simulated galaxies. Among the existing cosmological hydrodynamical simulations for the formation and evolution of galaxies, the IllustrisTNG simulations 2 have been particularly successful in reproducing a broad range of observational findings 34 . These include the galaxy color, stellar age, and metallicity trends at z ∼ 0 as function of galaxy stellar mass from SDSS as well as the qualitative characteristics of the stellar orbit distributions from the CALIFA Survey 35 . Among the flagship runs of the IllustrisTNG suite, the so-called TNG50 run 3, 4 is uniquely well-suited for (a) (b) Figure 5: Selection of simulated galaxies from the TNG50 simulation. Panel (a): Total stellar mass M * versus size R e of the simulated galaxies with the black selection box. Panel (b): For each of the 50 selected simulated galaxies, the stellar mass in the disk versus inner stellar halo component within 2 R e radius. The four components follow from the same division in phase-space of radius r versus circularity λ z as applied to NGC 1380. In both panels, the star symbol indicates the corresponding quantities for NGC 1380, while the cross symbols are five of the selected simulated galaxies located close to NGC 1380 in panel (b): they will serve as case studies.
our study, as (1) it reaches the numerical resolution of typical zoom-in simulations (with stellar particle of 8.5 × 10 4 M , so that the inner structural details of galaxies like thin disk, thick disk, bulge, and halo are well resolved, and (2) it still offers a large sample of galaxies (∼ 800 galaxies with M * > 10 10 M at z = 0) in a representative cosmological volume of about (50 Mpc) 3 that encompasses diverse environments.
In these simulations, halos and subhalos, and thus the basic properties of galaxies, are identified using the Friends-of-Friends (FoF) and SUBFIND algorithms 36, 37 . The past histories of all galaxies in the simulation can be followed via the merger trees, as constructed by the SUB-LINK code based on the baryonic mass of subhaloes 9 . When galaxies merge, the mass ratio of the galaxy-galaxy merger can be used as summary statistics: throughout, we will define such a ratio based on the stellar masses of the merging objects when the secondary progenitor reaches its maximum stellar mass. We selected galaxies from TNG50 at the z = 0 snapshot with total stellar masses and sizes similar to those of NGC 1380: 10 11 < M * < 4 × 10 11 M and 3 < R e /kpc < 8. This selection, indicated by a black box in Figure 5, results in 51 simulated galaxies, of which we excluded TNG50 411449 which has an on-going merger at z = 0. For each of the remaining 50 simulated galaxies, we measured their orbit-based disk, bulge, warm component, and inner stellar halo mass fractions within a 2 R e radius with their particle information from the simulations directly, adopting the same four-component division in r versus λ z as for NGC 1380. The five simulated galaxies having similar inner stellar halo mass as NGC 1380 will serve as cases studies. Different galaxies are identified by their subfind IDs in the TNG50 simulations as labeled in the figures.
Orbit integration of TNG50 galaxies. In the population orbit-superposition model 1 , we characterize the stellar orbits with two parameters, the radius r and circularity λ z : both are calculated as the averages of the corresponding quantities for particles stored with equal time-step along the orbits. In principle, for the simulated galaxies, we could directly take the instantaneous values of the stellar particles at z = 0. For a more fair and consistent comparison, we instead proceeded with an orbital integration. Namely, we freeze the stellar particles of a simulation at z = 0, construct its simulated potential by adding up stars, dark matter and gas mass distribution, and integrate 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.
We integrated the stellar orbits for the five case study galaxies, so that throughout the paper their stellar orbit properties are time-averaged values. However, the orbital integration is somewhat expensive to run for all the simulated selected galaxies. For the rest of the TNG50 sample, we therefore adopted an approximate method 38 , whereby averaging was 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 and it is good enough for our four-component separation. For the five case study galaxies, we have obtained the four component separation based on both the time and phase-space average methods, finding consistent results.
The separation into four stellar components of observed and simulated galaxies. We decomposed each galaxy into four components in the phase-space radius r versus λ z . The criteria of the separation were determined by considering the structures in the probability density p(r, λ z ), stellar age t(r, λ z ) and metallicity Z(t, λ z ) distributions, as well as the origin of these orbits as inferred from the simulations.
We first noticed that the dynamically cold orbits with λ z > 0.8 are usually younger and more metal-rich. In the simulated galaxies, these orbits represent a clean thin stellar disk. The dynamically hot orbits, with λ z < 0.5, usually show in both their probability density and metallicity distributions a transition from a central peak toward a plateau, as illustrated in panel (a) of Figure 6 for a typical simulated galaxy from the selected TNG50 sample, and panel (b) for NGC 1380 revealing a similar sharp drop plus plateau 3 Figure 6: Criteria for the four-component decomposition. Panels (a): For one of the simulated galaxies that is representative of the selected TNG50 sample, the probability density distribution and average metallicity, respectively, as function of radius r by summing all the orbits with λ z < 0.5 are shown. The vertical black dashed lines represent the effective radius R e . The black solid line indicates the radius r cut that separates the bulge from the inner stellar halo and is chosen where both the probability density and metallicity profiles flatten. Panel (b): The same as in panels (a) but for NGC 1380, with less smooth profiles due to model degeneracies, but still a visible flattening around r cut 3.6 kpc; the vertical red lines mark the extension of kinematic data r max . Panel (c): For all 50 simulated galaxies of the selected TNG50 sample, the distribution of r cut in units of kpc. Panel (d): For all the 50 simulated galaxies, the Pearson correlation coefficient between the stellar mass of the most massive merged satellite M * ,sat1 and that of the inner stellar halo M * ,halo(r<2Re) selected in phase-space for radius r cut < r < 2 R e and circularity λ z < λ z,cut . We adopt λ z,cut = 0.5 leading to the strongest correlation. Panel (e): similar to panel (d) but taking r cut = 3 kpc for all galaxies.
r cut in all the selected 50 simulated galaxies, we see from panel (c) that r cut is distributed around 3 kpc. We included all the orbits with r < r cut as the bulge component, except those on cold disk orbits. There are usually only few stars on cold disk orbits within r cut .
After separating the disk and bulge, the remaining orbits with r > r cut and λ z < 0.8 in the simulated galaxies seem mainly induced by mergers, even though stars with circularity closer to λ z = 0.8 can in principle be heated from the disk. We hence introduce a λ z,cut below which orbits with r > r cut are assigned to a so-called inner stellar halo component which is most strongly correlated with the most massive merger. To this end, we extract from all the 50 simulated galaxies the stellar mass of the most massive merged satellite M * ,sat1 as well as the mass of the inner stellar halo M * ,halo for different values of λ z,cut . Panel (d) of Figure 6 shows that the corresponding correlation is strongest for λ z,cut = 0.5, and similarly in panel (e) by taking r cut = 3 kpc for all galaxies. This is consistent with the scenario we learn from cases studies that the inner stellar halo is mainly induced by the most massive merger: the heated stars from the progenitor disk/bulge, deposited satellite stars, or stars formed during the massive merger(s) typically end up on orbits with r > r cut and λ z < 0.5.
While the latter defines the inner stellar halo component, the remaining orbits with r > r cut and 0.5 < λ z < 0.8 are assigned to a so-called warm component which can have complex and mixed origins. As seen in the galaxy simulations, a fraction of these stars in the warm component can also be induced by a massive merger: they can come from the progenitor disk being partially disrupted, or formed during the merger, but are rarely accreted from a massive merger as those stars typically end up on hot orbits. On the other hand, they could have been born warmer from more turbulent gas at higher redshift or heated by and accreted from minor mergers 39 .
Recovery of the four components. We have previously tested the ability of our population-orbit method to recover the age and metallicity properties of different orbital components, but so far only for simulated spiral galaxies without a significant inner stellar halo within the data coverage 1 .
Here we create mock data for TNG50 468590, which has similar inner stellar halo mass and disk mass as NGC 1380. The mock data are created with the same inclination angle ϑ = 76 • , similar data coverage and data quality as for NGC 1380.
We show in Figure 7 that our population-orbit method is able to fit the mock data in detail. Next, in panel (a) of Figure 8, we confirm that the probability density distribution p(λ z , r), age distribution t(λ z , r), and metallicity distribution Z(λ z , r) directly measured from the simulations are well matched by those extracted from the best-fitting population-orbit models. The only discrepancy is in the inner parts, where in the simulations the bulge is dominated by box orbits with λ z = 0 whereas in the model part of the box orbits are represented by tube orbits with non-zero in mainly the box orbits as the kinematic data do not extend far enough along the minor axis.  Figure 7: Population-orbit model to TNG50 468590. The first column shows the mock data from the simulated galaxy by projecting with similar inclination angle as NGC 1380, the second column the best-fit population-orbit model, and the third column the residuals of data minus model and divided by the errors assigned assuming similar quality data as for NGC 1380. From top to bottom: surface brightness, mean velocity, velocity dispersion, Gauss-Hermite coefficients h 3 and h 4 , as well as age and metallicity maps.   Figure 8: Verification of population-orbit method with TNG50 468590. Panel (a): From left to right, distributions in phase-space of radius r versus circularity λ z of orbital probability density p(r, λ z ), age t(r, λ z ), and metallicity Z(λ z , r), all are shown only within 2 R e . The top row shows the distributions directly from the simulations, which are well matched by those extracted from the best-fitting population-orbit models in the bottom row. Note that the true probability distribution P (r, λ z ) shown here is the same as what we shown in Figure 3, but with slightly different spatial coverage and unit, here 140 arcsec = 17 kpc with the mock galaxy placed at a distance of d = 25

20
Mpc. The dashed lines indicate our orbital-based division into four components: disk, bulge, warm component, and inner stellar halo. The corresponding mass fractions within a 2 R e radius are respectively 0.21, 0.43, 0.04, 0.32 from the simulations, which are well matched by 0.21, 0.40, 0.07, and 0.32 from the models. Panel (b): From top to bottom, stellar surface density, age, and metallicity profiles of the four components along the major axis. The left column shows the profiles directly coming from the simulations, which are well reproduced from those reconstructed from the models fitted to the mock data and shown in the right column.
but opposite values of λ z . This degeneracy might at least partially be the result of insufficient kinematic constraints along the minor axis. The latter, however, does not affect our division into four components: disk (λ z > 0.8), bulge (λ z < 0.8, r < r cut ), warm component (0.5 < λ z < 0.8, r > r cut ), and inner stellar halo (λ z < 0.5, r > r cut ). For TNG50 468590, the mass fractions within a 2 R e radius of the disk, bulge, warm component, and inner stellar halo are respectively 0.21, 0.43, 0.04, 0.32 from the simulations, which are well matched by 0.21, 0.40, 0.07, and 0.32 from its modeling.
For each of these components, we can then infer the stellar surface density, age and metallicity profiles along the major axis. Panels (b) of Figure 8 show that the profiles directly coming from the simulations are mostly reproduced from those reconstructed from the models fitted to the mock data.
Uncertainties on the mass of NGC 1380's inner stellar halo. We have shown that the mass fraction of the inner stellar halo is well recovered by our population-orbit model fitted to the mock data of a NGC 1380-like simulated galaxy. Whereas quantities in the simulations are massweighted, those for NGC 1380 are all r-band luminosity weighted.
The luminosity fractions for NGC 1380 of disk, bulge, warm and inner stellar halo components within 2 R e from our best-fit model are 0.19, 0.38, 0.17, and 0.26, respectively. Under the simplest assumption of our best-fit constant stellar mass-to-light ratio M * /L r = 2.7 ± 0.3 M /L from our dynamical models, we arrive at M * ,halo(r<2Re) = (3.4 ± 0.3) × 10 10 M , which we use throughout to infer the masses of the merged satellites. The uncertainty of M * /L r come from the 1σ scatter of the dynamical models.
In NGC 1380, M/L r varies from 4.8 ± 0.7 M /L in the bulge-dominated inner regions to 3.5±0.5 M /L in the inner-stellar-halo-dominated outer regions from stellar population synthesis by considering a variable IMF 40 , here the error is the 1σ scatter of M/L r in the corresponding regions. We have tested in previous work 7 that assuming constant mass-to-light ratio (even if it actually varies) to infer the gravitational potential does not affect the stellar orbit distribution in the dynamical model 7,41 , so that the luminosity fraction of the inner stellar halo is well represented by the adopted value of 0.26.     Case studies: inner stellar halo formation in TNG50 galaxies. We show here the formation of the inner stellar halos of other four simulated galaxies from TNG50, to support the findings obtained by studying TNG50 468590: these are TNG50 372754, 375073, 439099 and 117253. Like TNG50 468590, TNG50 372754 and TNG50 375073 have similar inner stellar halo mass as NGC 1380 (see Figure 5), but the three of them have different disk masses. TNG50 117253 and TNG50 439099 have less massive inner stellar halos.
We recall that we define the starting of a merger as the time when the centres of the galaxies get as close as 200 kpc for the first time; a merger is considered to be over by the time when the centers of galaxies have vanishing distance and do not depart from each other anymore.
TNG50 372754 experienced a massive merger with stellar mass ratio of ∼ 2 : 1 started at redshift z ∼ 0.4. Its evolution history is visualized in the top panels of Figure 9. Before the merger, the main progenitor was dominated by a cold disk and a compact bulge with a small fraction of inner stellar halo. During the massive merger, the disk of the main progenitor was heated to be on orbits of warm or inner stellar halo components. The satellite stars were accreted partly to the bulge, but mostly in the inner stellar halo. During the merger, stars were formed on the disk and partially heated to be warm. After the merger, newly-formed stars are on dynamically cold disk-like orbits, and persist as such until z = 0 because of the quiescent merger history after this event. The inner stellar halo mass of TNG50 372754 at z = 0 is 3.3 × 10 10 M , similar to that of NGC 1380, and ∼ 30 per cent of it was already there in the main progenitor before the merger event started at z = 0.4. ∼ 20 per cent was redistributed into the halo regions from the bulge and disk of the main progenitor through this merger and ∼ 50 per cent was accreted from the satellite galaxy. The contribution of other sources are negligible. In total 70 per cent of the inner stellar halo mass of TNG50 372754 was induced by the most massive merger.
TNG50 375073 had a massive merger with stellar mass ratio of ∼ 1 : 0.5 started at redshift z ∼ 1.3 as shown in middle panels of Figure 9. Its inner stellar halo was formed in a way similar to the two cases already described. It has an inner stellar halo mass of 3.2 × 10 10 M at z = 0: ∼ 14 per cent of it was already there in the main progenitor at z = 1.3; ∼ 38 per cent was redistributed into the halo regions from the bulge and disk of the main progenitor through this merger; ∼ 34 per cent was accreted from the satellite galaxy, and ∼ 12 per cent was formed during the merger event. In total 84 per cent of the inner stellar halo mass of TNG50 375073 was induced by the most massive merger.
TNG50 117253 accreted three massive satellites continuously from z = 2.5 to z = 1 as shown in lower panels of Figure 9. By considering a snapshot of this galaxy at z = 1.8, when the first merger ended, we can see that 12 per cent of the inner stellar halo mass was already in place at that time. Then the second and third 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 per cent of the inner stellar halo stars were accreted from the most massive satellite; ∼ 15 per cent were accreted from the second massive satellite, and ∼ 40 per cent were formed and induced during the mergers. At the end, TNG50 117253 has a total inner stellar halo mass of 2.2 × 10 10 M . The fraction produced by the most massive merger is less than 60 per cent.
TNG 439099 accreted two massive satellites: the most massive one was accreted at 1.5 < z < 1 together with a small satellite at almost the same time; later at 0.7 < z < 0.3 it accreted a slightly less massive one. We show the structure formation through these two merger events in Figure 10. TNG50 439099 has an inner stellar halo mass of 2.3 × 10 10 M at z = 0: ∼ 9 per cent of it was already in the main progenitor at z = 1.5; ∼ 75 per cent was formed during the first merger(s); ∼ 16 per cent of the inner stellar halo was induced during the second merger, mainly accreted from the satellite. The second merger heated the cold disk, which was reformed after the first merger ended, to be warm component.
The variety of the merger histories showcased by our five study cases may partially explain the scatter of the inner stellar halo mass vs. most massive satellite mass relation uncovered from TNG50 and shown in Figure 4. TNG50 468590 history is characterized by only one dominating merger that produced 90 per cent of the inner stellar halo: thus it is located above the median relation. In contrast, TNG50 117253 underwent multiple similarly massive mergers that all contributed to the inner stellar halo, and the most massive one only produced < 60 per cent of it: thus it is located below the median relation. TNG50 375073 and TNG50 439099 are on the average relation. TNG50 372754 had the most massive merger recently at 0.2 < z < 0.4, which heated the disk to be warm, but appeared not efficient in redistributing the disk stars into radial orbits (the inner stellar halo): this one is located above the average relation.
Note that in both TNG50 372754 and TNG50 439099, a massive merger occurring at low redshift efficiently heated the previous cold disk and produced a significant warm component: the stars in such warm components exhibit a wide range of stellar ages (Figure 11), and are on average younger than the stars on dynamical hot orbits with λ z ∼ 0. On the other hand, the warm component in NGC 1380 is as old as the stars on hot orbits, like those in TNG50 468590, 375073 and 117253. In all the five simulated galaxies we focused on, stars formed after the last massive mergers rarely ended on non-disk orbits. In all these cases, a sharp transition of stars from non-disk dominated orbits to disk dominated orbits can be clearly identified: therefore we can confidently associate this transition to the last massive merger ending time. In practice, the time when t| f disk =0.5 is a good proxy for the merger ending time, as shown in Figure 11.
Correlations between inner stellar halo of simulated galaxies and merged satellites. The five case studies show that the inner stellar halo is mainly induced by the most massive merger(s). This should be in general true for all galaxies as supported by the strong correlations we found between inner stellar halo mass and the most massive satellite(s) mass in panel (a) of Figure 4. If other processes, e.g., the combination of many small mergers, are also efficient in creating a massive  Figure 11: Timing of the last massive mergers of the five case-study simulated galaxies. The five panels show the five simulated galaxies selected for in-depth study. Similarly to panel (d) of Figure 4, in each panel, the bottom plots show the probability distribution of stellar orbits p(λ z , t) in the plane of circularity λ z versus stellar age t directly from the simulated galaxies. inner stellar halo, there should be galaxies located in the bottom-right corner of panel (a), thus diluting these strong correlations. And there are non of such points in our figure for the selected Illustris TNG50 sample.
The inner stellar halo in Figure 4 is separated by adopting λ z,cut = 0.5 and the standard transition radii r cut . The latter, however, is determined by eye and might introduce some arbitrariness. In Figure 12, we hence show the correlations from TNG50 by considering slightly different choices of r cut and λ z,cut for separating the inner stellar halo. The thick black lines are the fits to the correlations with the standard inner stellar halo, the same as shown in Figure 4. The thick red lines return the average relations adopting r cut = 3 kpc for all galaxies. The thick magenta lines use M * ,halo(r<2Re) + M * ,warm(r<2Re) in x-axis, which is equivalent to the case of taking λ z,cut = 0.8 for separating the inner stellar halo.
The correlations of M * ,halo(r<2Re) to the satellite(s) mass are almost identical irrespective of the choice for r cut . The vertical black and red lines mark the positions of NGC 1380 for the corresponding M * ,halo(r<2Re) . Note that by adopting r cut = 3 kpc instead of the transition radius of 3.6 kpc, M * ,halo(r<2Re) of NGC 1380 slightly increases. The satellite(s) mass for NGC 1380 we infer from these correlations are almost identical for these two options of r cut . The formation of these components are altered by the most massive mergers, but are unlikely to be mainly induced by them. The total accreted mass is, as expected, positively correlated with the total stellar mass of the galaxy 42, 43 . For comparison, we show in Figure 13 the correlations between the galaxy total stellar mass M * with the most massive satellite stellar mass in panel (a), with R = 0.54. The correlation with the mass of all satellites is shown in panel (c) and yields R = 0.68. These correlations are much weaker than the correlations between inner stellar halo mass M * ,halo(r<2Re) and the satellite(s) mass. At fixed stellar mass, the accreted stellar mass can vary significantly from galaxy to galaxy: on the other hand, the inner stellar halo mass can constrain the mass of accreted satellites also for single galaxies.
The merger in MGC 1380 is the most massive and oldest major merger discovered so far. In Figure 14, we contrast the most massive merger events in the only three galaxies that have been uncovered. The merger in the Andromeda galaxy happened only about 2 Gyr ago, and still substantial substructures, e.g, giant streams could be identified from the stellar halo. The most massive merger in the Milky Way is similarly old, but is minor with stellar mass of ∼ 4 × 10 8 M , which is only about two per cent of NGC 1380's progenitor satellite stellar mass.  Figure 14: NGC 1380's discovered major merger in comparison with past mergers in the Milky Way and the Andromeda galaxy. Gaia-Enceladus is thought to be the remains of a satellite with stellar mass ∼ 4 × 10 8 M , which merged with the Milky Way ∼ 10 Gyr ago 18, 19 . The giant stream in the inner stellar halo of the Andromeda galaxy (M31) indicates a major merger with a satellite of ∼ 2 × 10 10 M around ∼ 2 Gyr ago 21 . At a distance thirty times further than M31, our population-orbital modelling of the unresolved stars in NGC 1380, reveals that its inner stellar halo was induced by a major merger of ∼ 3 × 10 10 M that ended ∼ 10 Gyr ago.