2.1. Milky Way-like galaxies in TNG50

The IllustrisTNG simulations represent the next generation of cosmological hydrodynamical simulations designed to model galaxy formation and evolution within the framework of the $\Lambda$ CDM paradigm (Pillepich et al. Reference Pillepich2019; Nelson et al. Reference Nelson2019a). Building upon the foundations laid by the earlier Illustris simulations (Vogelsberger et al. Reference Vogelsberger2014a,Reference Vogelsbergerb; Genel et al. Reference Genel2014; Sijacki et al. Reference Sijacki2015), IllustrisTNG incorporates significant improvements, particularly in the modelling of AGN feedback, chemical enrichment, and the evolution of seed magnetic fields (Weinberger et al. Reference Weinberger2017; Pillepich et al. Reference Pillepich2018).

In this paper, we concentrate on the TNG50 simulation that operates within a periodic box with a size of L $_{\mathrm{box}}$ = 35 Mpc/h, containing 2 160 $^3$ gas elements and dark matter particles with mass resolutions of [0.85, 4.5] $\times 10^5 M_{\odot}$ , respectively. Moreover, TNG50 adopts the cosmological parameters specified by Planck Collaboration et al. (2016). This simulation provides a detailed and comprehensive framework for studying the formation and evolution of galaxies across cosmic time.

We examine a carefully selected sample of Milky Way-like galaxies from TNG50 identified and characterised in previous works (Emami et al. Reference Emami2020a,Reference Emamib; Emami et al. Reference Emami2022; Waters et al. Reference Waters2024). Our sample selection adheres to two primary criteria. First, we constrain the dark matter halo mass to fall within the range of (1–1.6) $\times 10^{12}\,\mathrm{M}_{\odot}$ , motivated by recent observational constraints on the dark matter halo mass of the Milky Way galaxy (Posti & Helmi Reference Posti and Helmi2019). Second, we limit our galaxy sample to rotationally supported galaxies identified in two steps, as the following.

As a first step, we compute the stellar net specific angular momentum vector, $\textbf{j}_{\textrm{net}}$ , for each galaxy in our sample:

(1) $\begin{equation}\textbf{j}_{\textrm{net}} \equiv \frac{\textbf{J}_{\textrm{tot}}}{M} = \frac{\sum_{i} m_i \textbf{r}_i \times \textbf{v}_i}{\sum_{i} m_i},\end{equation}$

where the index i refers to stellar particles. Aligning the z-axis with the direction of $\textbf{j}_{\textrm{net}}$ , we compute the inner product of each stellar particle’s angular momentum with the z-axis, given by $j_{z,i} = \textbf{j}_i \cdot \hat{\textbf{z}}$ . The orbital circularity parameter is then defined as

(2) $${\varepsilon _i} \equiv {{{j_{z,i}}} \over {{j_c}({E_i})}},\,\,\,\,{j_c}({E_i}) = {r_c}{v_c} = \sqrt {GM( \le {r_c}){r_c}} .$$

For each particle, we determine the radius of its corresponding circular orbit by equating its total energy to the specific energy of a circular orbit:

(3) $$E({r_c}) = {{GM( \le {r_c})} \over {2{r_c}}} + \phi ({r_c}),$$

where $ M(\!\leq r_c) $ denotes the mass enclosed within the circular orbit, and $ \phi(r_c) $ represents the gravitational potential at $ r_c $ , computed from the averaged radial profile of the total gravitational potential, accounting for contributions from stars, gas, dark matter, and the central black hole.

We define disk stars as those with $ \varepsilon_i \geq 0.7 $ . Furthermore, we restrict our sample to cases where at least 40% of the stars within a radial distance of 10 kpc from the centre exhibit disk-like properties. This criterion reduces our sample to 25 Milky Way-like galaxies.

2.2. H3 observation and selection functions

The H3 (Hectochelle in Halo at the High Resolution) Survey is an ongoing stellar spectroscopic survey that offers an unbiased measurement of stellar parameters (Conroy et al. Reference Conroy2019b, a). It provides spectro-photometric distances for approximately 2 $\times 10^5$ stars within the photometric magnitude range of $15 \lt r \lt 18$ , with a 3D heliocentric distance of d $_{\mathrm{helio}} \gt$ 3 kpc, $|b| \gt 40 ^{\circ}$ , and Dec $\gt -20^{\circ}$ . H3 survey outputs radial velocities, spectroscopic distances, [Fe/H], and [ $\alpha$ /Fe] abundances for the aforementioned stellar sample. By combining this data with the Gaia proper motion measurements, enables the determination of the full 6D phase-space information and 2D chemical-space information for these stars.

In a related study, Naidu et al. (Reference Naidu2020) focused on 5684 K giants from the H3 survey and conducted an extensive exploration of the structure of distant galaxies up to 50 kpc from the galactic centre. They analysed the scale-height dependence of the in situ stars. Naidu et al. (Reference Naidu2020) used a dedicated approach, combining the kinematics information as well as the metallicity and made a general class of in situ stars including the contributions from different parts such as the high- $\alpha$ disk, in situ halo, the metal-weak thick disk, Aleph (defined mainly based on their specific values of [Fe/H] $\gt$ −0.8 and $[\alpha$ /Fe] $\lt$ 0.27 together with some kinematic criteria such as the azimuthal component ( $-300\,\mathrm{km\,s}^{-1} \lt V_{\phi} \lt -175\,\mathrm{km\,s}^{-1} $ ) and the radial component ( $ \vert V_r \vert \lt 75\,\mathrm{km\,s}^{-1}$ ) of the stellar velocity), and unclassified disk debris.

Motivated by these observational findings, our study investigates in situ stars within a selected sample of Milky Way-like galaxies from the TNG50 simulation. We focus on assessing the scale height dependence of the in situ star fraction and contrasting our results with those from the H3 survey. While the H3 survey employed a combination of the kinematics and chemical cuts involving [Fe/H] and [ $\alpha$ /Fe] abundances, alongside photometric magnitude, in this first study, our approach to identifying in situ stars within the TNG50 simulations is based on the individual stellar birth locations being located within a given radial cut throughout the entire of the cosmic evolution. Following the notation of De Lucia & Blaizot (Reference De Lucia and Blaizot2007), the first progenitor is the one with the most massive assembly history of the main progenitor of the main subhalo in SubFind convention while the second progenitor is defined as the galaxy with the second massive history behind it. In this paper, we often refer to the first (f) and second progenitor (s). Furthermore, to facilitate comparison with the H3 findings, we add similar spatial selection functions and study the role of varying radial cut-offs to compile a sample of in situ stars in our galaxy dataset. Subsequently, we conduct a comparative analysis of the scale height profile of in situ stars with the H3 survey results, representing their distance from the galactic disk plane.

3.1. Inferring the in situ stars from TNG50 simulation

The in situ stellar component comprises stars born within the disk of their galactic host halo, as well as stars formed from streams of stripped gas originating from infalling satellites (Benson et al. Reference Benson, Lacey, Frenk, Baugh and Cole2004; Zolotov et al. Reference Zolotov2009; Purcell et al. Reference Purcell, Bullock and Kazantzidis2010; McCarthy et al. Reference McCarthy2012; Tissera et al. Reference Tissera, Beers, Carollo and Scannapieco2014; Pillepich et al. Reference Pillepich, Madau and Mayer2015; Rodriguez-Gomez et al. Reference Rodriguez-Gomez2016; Monachesi et al. Reference Monachesi2019). In situ stars born in the disk may be ejected to larger distances due various interactions. The amount of ejected mass may depend to some degree on subgrid physics and stellar feedback (see, for instance, Zolotov et al. Reference Zolotov2009; Cooper et al. Reference Cooper, Parry, Lowing, Cole and Frenk2015, and references therein). Consequently, the radial profile of the fraction of in situ stars varies among different simulations, ranging from negligible at larger distances ( $\gt5$ kpc) (Zolotov et al. Reference Zolotov2009; Pillepich et al. Reference Pillepich, Madau and Mayer2015) to dominant at larger radii (Font et al. Reference Font2011; Monachesi et al. Reference Monachesi2016b; Monachesi et al. Reference Monachesi2016a; Elias et al. Reference Elias2018; Monachesi et al. Reference Monachesi2019). On the contrary, accreted stars are born in satellite galaxies, with the satellite either bound to the main progenitor from the outset or accreted into the main progenitor at a later time (see, for example, Tissera et al. Reference Tissera, Beers, Carollo and Scannapieco2014; Monachesi et al. Reference Monachesi2019, and references therein).

In the following, we present our strategy for selecting the in situ stars for a sample of 25 MW-like galaxies in the TNG50 simulation. The in situ star fraction is defined as the ratio of in situ stars to the total number of stars at that location. For more details about our sample selection, please refer to Section 2.1.

We use merger trees constructed from the sublink algorithm (Rodriguez-Gomez et al. Reference Rodriguez-Gomez2015, Reference Rodriguez-Gomez2016) to identify the main progenitors of the subhalo, which are defined as the most massive progenitor of the target subhalo at any given cosmic time (De Lucia & Blaizot Reference De Lucia and Blaizot2007). Furthermore, we employ baryonic merger trees, tracking only the stellar components in our analysis. As previously mentioned, we exclusively consider stars belonging to the main progenitor. Specifically, we initiate with the stellar components bound to the main halo at zero redshift, $z = 0$ , and trace their birth locations backward in time. Subsequently, we ascertain whether each star was part of the main progenitor at its birth time. Finally, we determine the distance of the star from its host galaxy at its birth time and compare it with a radial threshold, denoted as $r_{\mathrm{cut}}$ . We employ four different selections for $r_{\mathrm{cut}}$ , encompassing either a constant value of 30 kpc comoving radius or three time-dependent distances defined as some fixed coefficient of the virial distance (hereafter $R_{\mathrm{vir}}$ ) at the star’s birth time; further details are provided below. If the star is situated within the specified threshold, we categorise it as an in situ star; otherwise, it is considered as being accreted onto the main halo.

3.2. Redshift evolution of the main progenitor

As previously mentioned, we use the SubFind main subhalos generated by the sublink algorithm to investigate the halo merger history for the first and second progenitors of all halos in our MW-like sample. Through this process, we determine the mass-ratio and merger time for each merger event. Following the standard methodology outlined in Rodriguez-Gomez et al. (Reference Rodriguez-Gomez2015, Reference Rodriguez-Gomez2016), the mass-ratio in a merger system is determined when the second progenitor reaches its maximum mass. However, estimating the precise merger time is more intricate, as the second progenitor may orbit around the first progenitor and subsequently lose a significant portion of its mass. In our analysis below, we compute all correlation coefficients at the point when the second progenitor achieves its maximum mass.

When analysing the influence of galaxy mergers on variables such as the star formation rate (SFR) and the comoving virial radius ( $R_{\mathrm{vir}}$ ), defined at the density threshold of 200 times the critical density, we categorise mergers into minor and major categories. Minor mergers are defined as those with a mass-ratio between 0.05 and 0.2, while major mergers are those with a mass-ratio above 0.2. However, for the computation of correlation coefficients, as described below, we focus exclusively on two scenarios: one involving mergers with a mass-ratio above 0.05 and the other involving mergers with a mass-ratio above 0.2.

3.3. The in situ stars from TNG50 vs H3

Having outlined our general methodology for inferring the in situ stars from the TNG50 simulation, we proceed to compare our theoretical predictions with the actual observational results from the H3 survey.

Figure 1 illustrates the $|Z|$ profile of the fraction of in situ stars in various halos within our sample. Hereafter Z denotes the vertical height away from the disk. To establish the Z direction, we perform a coordinate transformation from the simulation frame to a coordinate system wherein the Z direction aligns with the total angular momentum of the stellar disk. Once we define the Z-direction, we maintain flexibility in selecting the new $X-Y$ directions within the disk plane. Since the H3 survey provides results in the Heliocentric coordinate system, with the Sun located approximately 8 kpc away from the centre along the X direction, and given the absence of a real Sun in the TNG50 simulation, we repeat our analysis four times, placing the Sun in four orientations along the ( $\pm{\hat{i}}$ and $\pm{\hat{j}}$ ) directions. Subsequently, we average the final results obtained from the different orientations.

Figure 1. The |Z| profile depicts the fraction of in situ stars for different cut-offs, including (0.1, 0.2, 0.5) R $_{\mathrm{vir}}$ and 30 kpc, for a sample of 25 MW-like galaxies in the TNG50 simulation. In each panel, we have averaged over four different orientations for the Sun, located 8 kpc from the centre in the disk plane, spanning $\pm{\hat{i}}$ and $\pm{\hat{j}}$ . This averaging process helps reduce noise from various orientations. Overlaid on each panel, the solid green line represents the results from the H3 survey. Notably, there is a fair general agreement between the TNG50 results and the H3 survey, particularly for lower radial cuts. However, the level of agreement diminishes with increasing the threshold radius.

In each panel, solid lines of different colours [pink, blue, black] represent results using $r_{\mathrm{cut}} = [0.1, 0.2, 0.5] R_{\mathrm{vir}}$ , respectively, while the solid yellow line depicts results using $r_{\mathrm{cut}} = 30$ kpc. Additionally, the solid green line overlaid in each panel represents results from the H3 survey.

From the diagram, several observations can be made:

• • Increasing the radial cut-off does not drive the ratio of in situ stars to zero, indicating that in situ stars are present throughout most galaxies.

• • The $|Z|$ profile of the ratio of in situ stars varies among different galaxies. While some galaxies exhibit a significant decrease in $f_{\mathrm{in}}$ towards larger $|Z|$ s, others display a smoother profile.

• • Different galaxies exhibit varying sensitivity to different radial cut-offs. In some cases, increasing the radial cut-off notably enhances $f_{\mathrm{in}}$ , whereas in others, it demonstrates less sensitivity.

• • A constant cut-off of 30 kpc yields results similar to those obtained with a low radial cut of 0.1 $R_{\mathrm{vir}}$ . This observation is intriguing as the curvature of the profile remains largely consistent across most cases, despite the former choice being a constant independent of individual galaxy specifics or the exact birth time of stars.

• • In summary, a few galaxies [1, 11, 16, 18, 19, 23, 25] reasonably explain the trend in H3 results. Among these, only galaxies [1, 19, 25] match the observational results both in terms of the shape as well as the values of the in situ fraction regardless of the particular spatial cut for inferring the in situ fraction.

Figure 2 displays the Z profile of the fraction of in situ stars for all galaxies using two different radial cut-offs. The left panel corresponds to $r_{\mathrm{cut}} = 0.1 R_{\mathrm{vir}}$ , while the right panel depicts $r_{\mathrm{cut}} = 0.2 R_{\mathrm{vir}}$ . It is evident from the diagram that increasing $r_{\mathrm{cut}}$ leads to a greater deviation of the inferred in situ fraction from the observational results obtained from the H3 survey. Conversely, at $r_{\mathrm{cut}} = 0.1 R_{\mathrm{vir}}$ , most galaxies align well with the observational results. The lower cut-off appears more consistent with the constant cut-off used in previous literature (see for example Bonaca et al. Reference Bonaca, Conroy, Wetzel, Hopkins and Kereš2017).

Figure 2. The |Z $_{\mathrm{gal}}$ | profile illustrates the fraction of in situ stars to the total number of stars derived from the TNG50 simulation using two different radial cut-offs: $r{\mathrm{cut}} = 0.1 R_{\mathrm{vir}}$ (left) and $r_{\mathrm{cut}} = 0.2 R_{\mathrm{vir}}$ (right). Each panel encompasses the entirety of TNG50 results, with the observational results from the H3 survey overlaid on each diagram.

Several factors may contribute to the discrepancy in the in situ stellar fraction between H3 observations and TNG50 simulations. One possibility is numerical heating, which can artificially increase the scale height, leading to systematic deviations between theoretical predictions and observations. Alternatively, the observed differences may stem from variations in the merger histories of individual halos. Acknowledging these potential systematic effects as limitations of our comparison, we proceed in the following sections to investigate the relationship between the deviation in the in situ stellar fraction from the H3 results and the merger histories of the galaxies in our sample.

3.4. On the origin of in situ stars

There have been various hypotheses proposed regarding the origin of in situ stars in galaxies. Gao et al. (Reference Gao2010), Griffen et al. (Reference Griffen2018) suggested that the in situ component might originate from the collapse of cold gas, which is brought in from numerous gas-rich mergers. Conversely, they proposed that accreted stars could form from stellar sub-halos accreted onto the host galaxy. However, recent advancements, particularly through the analysis of Gaia DR1 and DR2, have significantly expanded our understanding beyond the solar neighbourhood, providing insights into the galactic halo as well. Gaia data suggests that the in situ component of the galactic halo might be associated with a population of heated disk stars (see for example Bonaca et al. Reference Bonaca, Conroy, Wetzel, Hopkins and Kereš2017; Haywood et al. Reference Haywood2018; Di Matteo et al. Reference Di Matteo2019; Belokurov et al. Reference Belokurov2020; Sestito et al. Reference Sestito2019, and references therein). These observational findings align well with N-body analyses (see for instance Zolotov et al. Reference Zolotov2010; Purcell et al. Reference Purcell, Bullock and Kazantzidis2010; Font et al. Reference Font2011; Qu et al. Reference Qu, Di Matteo, Lehnert, van Driel and Jog2011; McCarthy et al. Reference McCarthy2012; Jean-Baptiste et al. Reference Jean-Baptiste2017, and references therein).

The Toomre diagram is a powerful tool for distinguishing different stellar components of the Milky Way, particularly between the stellar disk and halo populations (see for example Bonaca et al. Reference Bonaca, Conroy, Wetzel, Hopkins and Kereš2017; Di Matteo et al. Reference Di Matteo2020, and references therein). Disk stars typically exhibit prograde motion with azimuthal velocities close to the Local Standard of Rest (LSR) (defined below), and their density gradually decreases as velocities deviate from the LSR, eventually reaching retrograde speeds. Since the relative number densities of disk and halo stars are shaped by past galaxy mergers, the Toomre diagram serves as an effective diagnostic for identifying merger remnants. When combined with metallicity and chemical abundance data, it provides a deeper link between stellar kinematics and population properties, offering insights into the formation history and evolutionary process of the Milky Way. Motivated by these insights, we use the Toomre diagram to investigate the redshift evolution of stars in various orbits, encompassing both members of the stellar halo and the stellar disk and investigate the role of galaxy mergers in altering the relative number of disk versus halo stellar components. As a part of our analysis, we construct the Toomre diagram for a subset of galaxies in our MW-like galaxy sample at redshift zero. Subsequently, we leverage this diagram across the entire redshift evolution to classify stars into halo-like and disk-like components. Furthermore, we explore the role of galaxy mergers in influencing the evolution of the halo-like and disk-like components.

Figure 3 illustrates the Toomre diagram at redshift zero for a subset of 4 galaxies in our sample, with all stars being included in the diagram. Each panel displays an overlaid dashed line representing the boundary between the stars in the disk-like and those on the halo-like orbits, defined as:

(4) $$\matrix{ {|V - {V_{{\rm{LSR}}}}| \ge {V_{{\rm{LSR}}}},\,\,\,\,\,\,{\rm{Halo}}} \cr {|V - {V_{{\rm{LSR}}}}| < {V_{{\rm{LSR}}}},\,\,\,\,\,\,{\rm{Disk}}} \cr } $$

where we define $V_{\mathrm{LSR}}$ as the local standard of rest (LSR) speed, which is equivalent to the circular speed, $V_{\mathrm{cir}}(r)$ , at the location of the Sun. Here, $V_{\mathrm{cir}}(r) \equiv \sqrt{G M_{\mathrm{tot}}(r)/r}$ , where $M_{\mathrm{tot}}(r)$ refers to the total mass interior to distance r including DM, star and gas particles. We consider the sun to be located at 8 kpc away from the galactic centre at redshift zero. The x-axis of the Toomre diagram depicts the azimuthal component of individual stars. It is inferred in two steps. First, we compute the total angular momentum of all of the stars. We subsequently infer the alignment of the angular momentum of each stellar particle with respect to this vector, as well as the component of their location along this direction. $V_{\phi,i}$ is then given as:

(5) $${V_{\phi ,i}} = {j_{z,i}}/\left( {{m_i}{R_{{\rm{perp}},i}}} \right),\,\,\,\,\,\,{R_{{\rm{perp}},i}} \equiv \sqrt {r_i^2 - z_i^2} $$

For convenience, we have removed the sub-index i in Figure 3. while the y-axis showcases the velocity orthogonal to the stellar disk angular momentum.

Figure 3. The Toomre diagram depicting the distribution of stars in four galaxies from our MW-like galaxy sample. The white contours in each panel (from the inner part to the outer part) refer to the top 10%, 30%, 50%, 70%, and 90% of stars, respectively. It is inferred that in each case there is a tail of stars on halo orbits, while the majority of stars are on the disk orbits. The dotted orange line in each panel shows the boundary between the halo and disk stars.

Overlaid in Figure 3, we have drawn contours representing the top 10%, 30%, 50%, 70%, and 90% of stars, ordered from the innermost to the outermost regions. It is evident that, across all galaxies, the top 90% contour reveals a tail of stars transitioning into halo orbits. However, the prominence of this tail at lower percentile levels varies depending on the individual galaxy. For instance, in galaxy 1, the Toomre diagram is highly compact, with nearly all stars—up to the 90% contour—remaining in the disk. In contrast, galaxy 14 exhibits the most extended stellar tail in halo orbits, including a population with retrograde motion. Galaxies 3 and 25 present intermediate cases, where a significant number of stars are distributed in both the disk and halo.

The selection of these galaxies is based on their in situ fraction profiles from Figure 1. Galaxies 1 and 25 exhibit profiles closely resembling the H3 observations, while galaxies 3 and 14 show progressively greater deviations from the observations as the spatial cut increases. Galaxy 1 exhibits a highly compact Toomre diagram, consistent with its strong similarity to the H3 survey and the presence of very few stars on halo orbits. In contrast, galaxy 14 displays the most extended Toomre diagram, with a prominent tail of stars in halo orbits. This trend aligns with the increased deviation between the in situ stellar fraction in TNG50 and the H3 survey at larger spatial cuts. Interestingly, the Toomre diagrams of galaxies 3 and 25 appear very similar, despite notable differences in their in situ stellar fraction trends. In the case of galaxy 3, deviations emerge at larger spatial cuts, whereas in galaxy 25, the agreement between theory and observation remains robust regardless of the spatial cut used to define in situ stars. This distinction is reasonable, as the Toomre diagram represents stellar distribution in velocity space, while the in situ stellar fraction is defined based on spatial coordinates. Given the possible limitation of connecting these two space, moving forward we will use alternative metrics to quantify the observed deviations between the MW-like galaxies from the TNG50 simulations and the H3 observation.

As $V_{\mathrm{LSR}}$ speed is essential in the distinction between the halo-like and disk-like stars, in Figure 4 we present the redshift evolution of this quantity, defined as the circular speed $V_{\mathrm{cir}}$ evaluated at $r = A_{0} R_{\mathrm{h}}(z)$ . Here, $R_{\mathrm{h}}(z)$ represents the redshift-dependent half-light radius, and $A_{0}$ is a constant determined such that $r = 8$ kpc at redshift zero. Thus, we have $A_{0} = 8/R_{\mathrm{h}}(z=0)$ . The choice of 8 kpc ensures that the $V_{\mathrm{LSR}}$ is computed at the Sun’s location at present.

Figure 4. The redshift evolution of the local standard of rest velocity for 25 MW-like galaxies in our sample. The velocity is in the unit of km/s.

The key takeaway from Figure 4 is the universality of the $V_{\mathrm{LSR}}$ as a function of redshift for different galaxies in our galaxy sample.

Figure 5 illustrates the redshift evolution of the ratio of halo (denoted by diamond-red markers) and disk (represented by grey-circle markers) in situ stars as identified at redshift zero defined using spatial cut at 0.2 $R_{\mathrm{vir}}$ . It is important to note that the choice of 0.2 $R_{\mathrm{vir}}$ is motivated by the need to include a larger number of stars while effectively capturing the impact of galaxy mergers on the orbits of both disk and halo stars. While the in situ stars are defined using a spatial cut at 0.2 $R_{\mathrm{vir}}$ , we only show the evolutionary trajectory of their number at r = 17–23 kpc at zero redshifts, to the total number of in situ stars at redshift zero within the same radial range. We have chosen this range, r = 17–23 kpc, as according to Figure 1, the in situ stellar fraction diminishes significantly in this interval. Furthermore, to simplify delivering the picture, we only present the results for a case with the Sun located at 8 kpc in the $+X$ direction.

Figure 5. The redshift evolution of the fraction of the number of disk (represented by grey-circle) and halo (depicted by red-diamond) in situ stars between r = 17–23 kpc to the total number of in situ stars at redshift zero with a cut-off of 0.2 $R_{\mathrm{vir}}$ . Additionally, minor and major mergers are indicated by dashed-plum and solid-purple lines, respectively. The thresholds for minor and major mergers are set as 0.05–0.2 and above 0.2, respectively. It is seen that both the major and minor mergers change the slope of the evolution of the fraction of in situ stars.

For each galaxy, we begin with the selected in situ stars at zero redshift, identified based on the criteria outlined in Section 3.3 and trace them back in time. At each redshift, we compute the $V_{\mathrm{LSR}}$ and categorise the stars into halo (diamond-red) and disk (grey-circle) components. Additionally, major mergers are depicted by solid-purple lines, while minor mergers are represented by dashed-plum lines.

The two ratios are complementary as expected. It is evident that stars on disk orbits begin with a negligible fraction at high redshifts and gradually increase in their percentage towards lower redshifts. Ultimately, at redshift zero, the in situ stellar fraction is predominantly dominated by stars in disk-like orbits, which aligns with the predominant stellar disk characteristic of MW-like galaxies. Galaxy mergers play a role in altering the growth/decline slope for both halo and disk stars.

Figure 6 illustrates the redshift evolution of the in situ stellar mass in disk orbits (left panel) and halo orbits (right panel). The original in situ stars are traced back in time, and in each snapshot, the stellar mass in the disk and halo orbits is inferred. It is observed that throughout the galaxy’s evolution, the total stellar mass in both disk and halo orbits increases, with the final mass in disk stars being dominant over their values in halo orbits.

Figure 6. The redshift evolution of the total mass of the disk in situ stars (left panel) and halo in situ stars (right panel), defined with a spatial cut-off of 0.2 $R_{\mathrm{vir}}$ .

3.5. Redshift Evolution of R $_{\mathrm{vir}}$ and SFR

With the merger history of the first progenitor elucidated, we now delve into the redshift evolution of two key parameters in our analysis: the comoving virial radius $R_{\mathrm{vir}}$ and the star formation rate (SFR). The comoving virial radius is pertinent as it influences the radial cut-off used to infer the in situ stars at each redshift. On the other hand, the SFR is linked to the in situ stellar budget at any given time. Both of these quantities are impacted by merger events.

Figure 7 presents the redshift evolution of $R_{\mathrm{vir}}$ and the SFR for our entire MW-like sample. In each diagram, the time of major mergers is indicated by a solid purple line, while minor mergers are represented by dashed plum lines.

Figure 7. Redshift evolution of the comoving virial radius (R ${\mathrm{vir}}$ ) (left) and the star formation rate (SFR) (right) for our sample of MW-like galaxies. Vertical solid lines (purple) indicate major mergers, while dashed lines (plum) represent minor mergers. It is observed that both R ${\mathrm{vir}}$ and SFR are enhanced during merger events.

Throughout the galaxy’s evolutionary trajectory, instances arise where either or both $R_{\mathrm{vir}}$ and SFR change their amplitude and slope. While there are various reasons for these slope variations, e.g. owing to galaxy-galaxy mergers, due to pseudo-evolution owing to cosmic evolution, below we mainly focus on their impact on the evolution of in situ star fraction. More specifically, we embark on an analysis to explore the correlation between the evolution of these quantities and the deviation between the in situ fraction inferred from the TNG50 simulations and the H3 survey.

3.5.1. Quantitative comparison between TNG50 and H3 results

Here we construct a more quantitative metric to facilitate the comparison between the scale-height distribution of the in situ stars from the H3 survey and the TNG50. It is noted from each panel in Figure 1 that the in situ fraction from H3 survey has two important features. Its sharp slope at smaller $\vert Z \vert$ , i.e. $\vert Z \vert:$ (5–7.5) kpc, followed by its amplitude decline at larger distances, i.e. $\vert Z \vert:$ (17–20) kpc. To make an in-depth comparison between the theory and observation we include both of the amplitude and the slope (i.e. the derivative of the in situ star fraction) inside our metric. We have made a dedicated test analysis of either dropping the amplitude difference or the slope difference from the metric and realised that in both cases the strength of the subsequent correlation functions (as will be extensively discussed in Section 4) with both of the internal and external drivers diminishes. Having this pointed out, we suggest the following quantitative metric for comparing the H3 survey with TNG50:

(6) $\begin{align}\Delta f^{j}_{\mathrm{eff}} \equiv & \bigg\vert \frac{ \sum _{i = 1}^{\mathrm{N_1}} \big{(} f^{j}_{\mathrm{in}}(i) - f_{\mathrm{H3}}(i) \big{)} }{2\mathrm{N_1}} \bigg\vert + \bigg\vert \frac{\sum _{i = 1}^{\mathrm{N_2}} \big{(} f'^{j}_{\mathrm{in}}(i) - f'_{\mathrm{H3}}(i)\big{)}}{\mathrm{N_2}} \bigg\vert \nonumber\\\equiv & \Delta f^{j}_{\mathrm{Amp}} + \Delta f^{j}_{\mathrm{Slope}},\quad j = (0.1, 0.2, 0.5).\end{align}$

here $f^{j}_{\mathrm{in}}(i)$ refers to the amplitude of the fraction of the in situ stars from the TNG50. The upper index describes the spatial cut in inferring the in situ stars from simulations, while the sub-index i presents the scale height that we have evaluated the in situ star fraction as outlined below. $f_{\mathrm{H3}}(i)$ describes the fraction of the in situ stars from the H3 survey located in ith bin. We split the scale height between $\vert Z \vert:$ (17–20) kpc to $N_1 = 23$ linear bins and evaluate the difference between the TNG50 and the H3 in each of these points. We then compute the mean value of the amplitude difference including all of these $N_1$ points.

In the second term we infer the difference between the derivative of the in situ fraction with respect to the scale height. Here we cover the scale height in the range 5–7.5 kpc and split this range to $N_2$ = 10 points. $f'^{j}_{\mathrm{in}}(i)$ refers to the derivative of the fraction of in situ star defined with the spatial cut at j with respect to the scale height, while $f'_{\mathrm{H3}}(i)$ describes the derivative of the in situ fraction from H3 survey evaluated at the bin ith.

Since the theoretically inferred in situ fraction depends on the actual cut-off, we may want to try different deviations at the level of 0.1 $R_{\textrm{vir}}$ , 0.2 $R_{\textrm{vir}}$ , and 0.5 $R_{\textrm{vir}}$ . Using $\Delta f^{j}_{\mathrm{eff}}$ , we can associate one number to every galaxy and use it as a proxy to judge how close the simulated galaxy is to the observations. Another advantage of this metric is that we can make an automated way to accurately infer the proximity between the theory and H3 observation.

Below, we extensively make use of $\Delta f^j_{\mathrm{eff}}$ when we connect the above deviations to the merger history and the actual properties of individual halos. Our goal is to find some possible correlations between how small/large is $\Delta f^j_{\mathrm{eff}}$ and the closeness of different MW-like galaxies chosen from the TNG and the actual MW from the H3.

4.1. Internal drivers for the in situ fraction

As discussed in Section 3.5, $R_{\mathrm{vir}}$ and the SFR are fundamental internal parameters with direct implications for galaxy evolution. Building upon this insight, we investigate the correlation matrix between these parameters and the deviation between the in situ fraction from the TNG50 simulation and H3 survey as given in Equation (6). Given the evolutionary nature of these internal drivers over a galaxy’s lifespan, we segment the data into redshift bins to analyse their significance at different epochs. Specifically, we compute the mean, median, and slope of these parameters using two distinct approaches: either from an initial redshift to $z_{\mathrm{cut}}$ (referred to as H), or from $z_{\mathrm{cut}}$ to redshift zero (named as L). Subsequently, we assess the correlation between these metrics and $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ derived from Equation (6).

We conduct an extensive exploratory analysis to assess the significance of different combinations and operations on $R_{\mathrm{vir}}$ and SFR in our study, constructing correlation matrices using both the mean and median of these parameters, as well as their derivatives. We also incorporate both the original and the absolute-valued derivatives.

Within each correlation matrix, the relative importance of different metrics fluctuates across $z_{\mathrm{cut}}$ values, highlighting the role of both amplitude and slope in capturing deviations between theory and observation. Notably, the mean SFR exhibits a strong correlation, while the slope of $R_{\mathrm{vir}}$ proves to be more significant than its value. Furthermore, the median slope of both SFR and $R_{\mathrm{vir}}$ emerges as more influential than their mean counterparts, suggesting that gradual variations in these quantities are more impactful than outliers.

Interestingly, for $\frac{dR_{\mathrm{vir}}}{dz}$ , the absolute-valued median exhibits stronger correlations, whereas for SFR and $\frac{d \mathrm{SFR}}{dz} \Big|_H$ , the mean and median absolute values show weaker correlations. This finding led us to adopt a flexible approach, selecting either the mean or median, as well as either the parameter value or its slope in our analysis.

Having this pointed out, we keep the following discussions concise and present only the most significant parameter combinations in our correlation matrix. This includes: $\widetilde{\frac{dR_{\mathrm{vir}}}{dz}}|_H$ , $\widetilde{\frac{dR_{\mathrm{vir}}}{dz}}|_L$ , $\overline{SFR}_L$ , $\widetilde{\frac{d SFR}{dz}}|_H$ , and $\widetilde{\frac{d SFR}{dz}}|_L$ , The median and mean of the X quantity are denoted by $\tilde{\mathrm{X}}$ and $\overline{\mathrm{X}}$ , respectively.

Figure 8 illustrates the correlation coefficient matrix for both the H and L redshift segments across varying $z_{\mathrm{cut}}$ values. Each row represents a spatial cut, progressively increasing from $0.1 R_{\mathrm{vir}}$ to $0.5 R_{\mathrm{vir}}$ . While there are discernible variations in correlation coefficients between the smallest ( $0.1 R_{\mathrm{vir}}$ ) and largest ( $0.5 R_{\mathrm{vir}}$ ) cuts, overall trends, including magnitude and sign, remain consistent.

Figure 8. The correlation coefficient of $\widetilde{|\frac{dR_{\mathrm{vir}}}{dz}|}_H$ , $\widetilde{|\frac{dR_{\mathrm{vir}}}{dz}|}_L$ , $\overline{\mathrm{SFR}}_L$ , $\widetilde{\frac{d \mathrm{SFR}}{dz}}|_H$ and $\widetilde{|\frac{d \mathrm{SFR}}{dz}|}_L$ with $\Delta f^j_{\mathrm{eff}}$ . Different matrices from top to bottom correspond to various spatial cuts of $j = (0.1, 0.2, 0.5)$ , respectively.

It is evident that in most cases, $\widetilde{|\frac{dR_{\mathrm{vir}}}{dz}|}_H$ is positively correlated with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ . To elucidate the underlying reason for this positive correlation, we delve deeper into Figure 7, where it becomes apparent that $R_{\mathrm{vir}}$ undergoes two distinct evolutionary phases: a growth phase followed by a contracting phase. The transition point from the former to the latter phase varies across different galaxies. An earlier rapid growth expands the boundaries of $R_{\mathrm{vir}}$ to larger distances, consequently increasing the fraction of in situ stars and thereby elevating the value of $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ . Moreover, Figure 8 suggests that for a very early redshift cut, i.e., the higher values of $z_{\mathrm{cut}}$ , larger spatial cuts are more correlated than the lower values. We posit that this phenomenon is associated with the inclusion of the transitional phase from an increasing to diminishing evolutionary phase of $R_{\mathrm{vir}}$ .

At lower values of $z_{\mathrm{cut}}$ , $\widetilde{|\frac{dR_{\mathrm{vir}}}{dz}|}_L$ is anti-correlated with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ , indicating that as $R_{\mathrm{vir}}$ diminishes, the spatial boundaries contract, resulting in a more concentrated distribution of stars and thus lower values of $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ . Conversely, at higher values of $z_{\mathrm{cut}}$ , this anti-correlation transitions to a positive correlation. We interpret this shift as a consequence of changes in the slope of $R_{\mathrm{vir}}$ corresponding to increasing $z_{\mathrm{cut}}$ .

Figure 8 illustrates a positive correlation between $\overline{\mathrm{SFR}}_L$ and $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ . This correlation signifies the potential for higher $\overline{\mathrm{SFR}}_L$ values to stimulate star formation across various regions of galaxies. Notably, the strength of this correlation fluctuates with $z_{\mathrm{cut}}$ : stronger correlations are evident at smaller spatial cuts for lower $z_{\mathrm{cut}}$ values, while stronger correlations occur at larger spatial cuts for higher $z_{\mathrm{cut}}$ values. We propose that this behaviour is driven by the evolutionary dynamics of galaxies. Due to the higher surface density and the Schmidt-Kennicutt relation, stars are more likely to be concentrated near the centre at low redshift, resulting in stronger correlations for smaller spatial cuts. In contrast, at higher redshifts, star formation occurs at larger distances, thereby enhancing correlations at larger spatial cuts.

Comparing the correlation coefficients of $ \widetilde{\left| \frac{dR_{\mathrm{vir}}}{dz} \right|}_H $ with $ \widetilde{\frac{d \mathrm{SFR}}{dz}}|_H $ , it is evident that they closely mirror each other in most cases. The key distinction lies in the former being the median of the absolute values, whereas the latter is simply the median. This difference leads to a reversal in the sign of the correlation coefficient: $ \widetilde{\left| \frac{dR_{\mathrm{vir}}}{dz} \right|}_H $ exhibits only positive correlations, while $ \widetilde{\frac{d \mathrm{SFR}}{dz}}|_H $ predominantly shows negative correlations. Despite this sign difference, their correlation magnitudes remain fairly similar.

In summary, the correlation coefficient structure exhibits remarkable consistency between $R_{\mathrm{vir}}$ and SFR. In most cases, the median of the absolute valued derivative shows a positive correlation with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ , underscoring the role of enhancements in these values in amplifying the in situ stellar fraction across galaxies.

Figure 9 showcases the cross-correlation between various internal metrics and $z_{\mathrm{cut}}$ . The top and bottom rows correspond to $\widetilde{\frac{dR_{\mathrm{vir}}}{dz}}|_H$ and $\widetilde{\frac{dR{\mathrm{vir}}}{dz}}|_L$ , respectively, in conjunction with stellar-driven metrics. Within each row, from left to right, we demonstrate the correlation of $\widetilde{|\frac{dR{\mathrm{vir}}}{dz}|}_{i = (H,L)}$ with $\overline{\mathrm{SFR}}_L$ , $\widetilde{\frac{d \mathrm{SFR}}{dz}}|_H$ , and $\widetilde{|\frac{d \mathrm{SFR}}{dz}|}_L$ , respectively. The presence of non-zero cross-correlations between different components suggests their shared contributions to $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ .

Figure 9. The cross-correlation coefficient between different internal metrics. In top/bottom row, from left to right, we present the correlation of $\widetilde{\frac{dR_{\mathrm{vir}}}{dz}}|_H$ / $\widetilde{\frac{dR_{\mathrm{vir}}}{dz}}|_L$ with $\overline{SFR}_L$ , $\widetilde{\frac{d SFR}{dz}}|_H$ and $\widetilde{\frac{d SFR}{dz}}|_L$ , respectively.

4.2. External drivers for the in situ fraction

During the evolutionary trajectory of galaxies, numerous external events, including mergers and disruptions, play pivotal roles in shaping their dynamics. These events significantly influence star formation rates and redistribute stars within galaxies. Our objective here is to investigate how galaxy mergers impact the distribution of in situ stars.

Figure 10 illustrates the influence of galaxy mergers on the creation and distribution of in situ stars in two galaxies (depicted in the top and bottom rows) of our galaxy sample. The left panels depict the birth location versus the current location, while the right panels display the look-back time (TLB) versus the birth location of the in situ stars located between 0.1–0.2 $R_{\mathrm{vir}}$ . The dotted-dashed-gold and dashed-white lines denote minor and major mergers, respectively. We have consistently computed the look-back time for all of the stars of interest as well as for mergers. Focusing on the right panels of Figure 10, it is evident that galaxy mergers contribute to the formation of high-density peaks in the stellar distribution, appearing as bright regions in the kernel density estimator (KDE). For instance, in the top (bottom) panel, a major merger occurring at $T_{\mathrm{LB}} \simeq $ 7.5 (11.5) Gyr results in the formation of a distinct density peak. Additionally, in the bottom panel, the KDE shape is noticeably altered following a minor merger at $T_{\mathrm{LB}} \simeq $ 6.5 Gyr, highlighting the impact of mergers on the stellar distribution.

Figure 10. Rows present the birth location vs the current location (left panels) as well as the look-back time ( $T_{\mathrm{LB}}$ ) vs the birth location (right panels) of the in situ stars located between the 0.1 $R_{\mathrm{vir}}$ and 0.2 $R_{\mathrm{vir}}$ in two galaxies of the sample analysed in this work. The dotted-dashed-gold and dashed-white lines refer to the minor and major mergers, respectively. Both minor and major mergers are important in elevating star formation and producing in situ stars. In addition, in situ stars get closer to the centre after their birth.

Both minor and major mergers play a significant role in enhancing star formation and generating in situ stars. Moreover, the analysis suggests that in situ stars tend to migrate closer to the centre following their birth. Motivated by this, in the following, we compile a comprehensive list of potentially influential parameters and analyse their correlation coefficients with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ . We provide detailed definitions of each parameter and elucidate their respective contributions to the observed deviations between the TNG50 and H3 datasets.

$\bullet$ Effective mass-ratio: One of the crucial parameters significantly influencing $\Delta f^j_{\mathrm{eff}}$ , $j = (0.1, 0.2, 0.5)$ , is the mass-ratio of galaxy mergers. Motivated by this insight, we introduce a metric termed the effective mass-ratio of galaxy mergers, $\overline{\mathrm{MR}}$ , which encompasses both stellar mass and star-forming gas mass-ratios, where we define star-forming gas as cells with instantaneous star formation, identified by selecting cells with a positive ‘StarFormationRate’ in their gas field. While the role of the stellar mass-ratio, between first and second progenitors, is relatively well-understood, the influence of star-forming gas warrants further investigation. We posit that a gas-rich galaxy merger collectively enhances star formation. Consequently, we incorporate this additional contribution into our metric and observe a consistent correlation coefficient across all cases.

For each of these contributions, we compute the average mass-ratio above a specified threshold, denoted as $f_{\mathrm{MM}}$ , after a designated redshift threshold. In subsequent analyses, we explore the impact of two values for the mass-ratio threshold, specifically $f_{\mathrm{MM}}=(0.05, 0.2)$ . Additionally, through exploratory investigations, we determine that a redshift threshold of $z \leq 4$ effectively provides a better correlation, as the majority of galaxy evolution occurs below this threshold. Furthermore, it is observed that this threshold is consistent with the mass-weighted mass-ratio as defined in Appendix A.

$\overline{\mathrm{MR}}$ is defined as:

(7) $\begin{equation}\overline{\mathrm{MR}} = \frac{1}{N}\sum_{i} ^{N} \left(\mathrm{MR}_{*} + \mathrm{MR}_{g} \right)_i,\end{equation}$

where the summation index i is over the total number of the mergers with the mass-ration above $f_{\mathrm{MM}}$ . In each merger event, we identify the first and the second progenitors with $M^f_i$ and $M^s_i$ , respectively. Furthermore, $\mathrm{MR}_{j} \equiv M^s_i/M^f_i$ , $j = *$ , and $j= g$ denote the mass-ratio of stars to star-forming gas, respectively.

$\bullet$ Mean distance: The second crucial parameter, correlated with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ , pertains to the impact parameter in galaxy mergers. In the event of a galaxy merger, the distance between the first and second galaxy significantly influences the redistribution of stars. Our analysis involves computing the average distance, hereafter referred to as $\overline{D}$ , between the first and second galaxy of any galaxy mergers with mass-ratio above $f_{\mathrm{MM}}$ happening below the $z \leq 4$ . The distance is inferred at the time when the second progenitor gets to the peak of its mass. Since we do not have enough time resolution in the cosmological simulations, it is common to approximate this snapshot, or its subsequent snapshot, at the time of the merger.

$\bullet$ Effective redshift: The third critical parameter influencing $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ , is the timing of galaxy mergers. The occurrence of a galaxy merger during the primary epoch of galaxy evolution may exert a more pronounced effect on the stellar profile compared to mergers happening later. We establish an effective redshift of mergers, denoted as $z_{\mathrm{merg}}$ , for mergers transpiring below $z \leq 4$ and surpassing a certain merger fraction threshold $f_{\mathrm{MM}}$ . This effective redshift is calculated as $z_{\mathrm{merg}} = (1/N)\sum_{i}^{N} z_i$ , where N represents the total number of galaxy mergers, while $z_i$ describes the redshift associated with the ith merger event.

$\bullet$ Number of mergers: The fourth pivotal parameter is the total number of mergers with mass-ratios exceeding $f_{\mathrm{MM}}$ occurring below $z \leq 4$ . We denote this quantity as $N_{\mathrm{merg}}$ and explore its correlation with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ .

$\bullet$ Mean spin ratio: The fifth pivotal parameter associated with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ , is the spin ratio between the first and second progenitors at the moment of a merger. Where the halo spin is defined as the total angular momentum vector, computed as the mass-weighted sum of the product of a particle’s coordinate and velocity for all particles or gas cells associated with the halo. From a physical perspective, in the event of a galaxy merger, both the spin magnitude and the spin alignment between the first and second progenitor play a crucial role. The spin magnitude, or more precisely the spin ratio, encodes how the stellar distribution is expected to evolve post-merger. A higher spin ratio between the second and first progenitor suggests that stars are more likely to be scattered farther after the merger, increasing the deviation between the in situ fraction in TNG50 and the H3 survey.

Spin alignment is also expected to be a key factor, particularly in distinguishing between face-on and edge-on merger configurations. We define the average value of this spin ratio, incorporating all galaxy mergers with mass-ratios exceeding $f_{\mathrm{MM}}$ occurring below $z \leq 4$ throughout the merger history of galaxies, and denote this quantity as $\overline{\left(\frac{S2}{S1}\right)}$ . Our investigation revealed that taking the logarithm of this quantity yields a stronger correlation. Therefore, we employ $\log{(\overline{\frac{S2}{S1}})}$ when computing the correlation coefficients.

$\bullet$ Maximum spin alignment: The final critical parameter associated with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ , is the spin alignment between the first and second progenitors. We differentiate between the impact of face-on and edge-on galaxy mergers on the redistribution of stars across galaxies. Our investigation has demonstrated that the maximum alignment, derived from a subset of galaxy mergers meeting the aforementioned criteria, yields a stronger correlation coefficient. Furthermore, it has been observed that the logarithm of this quantity provides an even higher correlation. We define this parameter as $\log{(S1 \cdot S2)_{\mathrm{M}}}$ .

4.2.1. Corr(External driver, $\Delta f^j_{\mathrm{eff}}$ )

Figure 11 illustrates the correlation function between individual external parameters and $\Delta f^j_{\mathrm{eff}}$ . The first and second rows present the correlation coefficients for $f_{\mathrm{MM}} = 0.05$ and $f_{\mathrm{MM}} = 0.2$ , respectively. Within each row, different colours represent distinct spatial cuts: green, blue, and pink correspond to 0.1 $R_{\mathrm{vir}}$ , 0.2 $R_{\mathrm{vir}}$ , and 0.5 $R_{\mathrm{vir}}$ , respectively.

Figure 11. The correlation coefficient between parameters associated with galaxy mergers and $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ . The top row presents the correlation coefficient for the galaxy mergers above the effective mass-ratio 0.05, while the bottom row depicts the correlation coefficient for mergers with mass-ratio above 0.2.

From the diagram, it is evident that all external drivers are positively correlated with $\Delta f^j_{\mathrm{eff}}$ , although the strength of correlation coefficients varies across different spatial cuts and mass-ratio merger thresholds. It is anticipated that higher values of $\overline{\mathrm{MR}}$ and $N_{\mathrm{merg}}$ result in larger $\Delta f^j_{\mathrm{eff}}$ . It is also seen that $\overline{D}$ positively correlates with $\Delta f^j_{\mathrm{eff}}$ . We argue that mergers with larger $\overline{D}$ , unless those with very close to equal mass ratios, will involve galaxies with multiple encounters with the first progenitor before the actual merger. This may trigger both star formation as well as further enhance the scattering in the stellar population. Furthermore, we speculate that galaxy mergers with larger impact parameters may also leave star-forming gas at larger distances, causing them to form stars at larger radii.

$z_{\mathrm{merg}}$ is also positively correlated with $\Delta f^j_{\mathrm{eff}}$ . Our analysis suggests that having galaxy mergers at higher redshifts, albeit below $z=4$ , is crucial for transporting star-forming gas to the outer regions of galaxies, thus facilitating star formation at subsequent redshifts.

The highest correlation coefficient at $f_{\mathrm{MM}} = 0.05$ is attributed to $N_{\mathrm{merg}}$ , suggesting its strongest association with $\Delta f^j_{\mathrm{eff}}$ . Nevertheless, both $\overline{\mathrm{MR}}$ and $z_{\mathrm{merg}}$ also exhibit significant importance.

At $f_{\mathrm{MM}} = 0.05$ , the lowest correlation coefficient is linked to the spin-driven quantities, although this changes at $f_{\mathrm{MM}} = 0.2$ . At this level, the correlation coefficients for all quantities at 0.2 $R_{\mathrm{vir}}$ increase, while some of them decrease for 0.1 $R_{\mathrm{vir}}$ . As shown in Figure 11, more aligned mergers, where the progenitor spins are nearly parallel (i.e., face-on mergers), exhibit a higher correlation coefficient with $\Delta f^j_{\mathrm{eff}}, j = (0.1, 0.2, 0.5)$ . A deeper investigation with higher time resolution would be valuable to quantitatively assess the role of face-on vs. edge-on mergers in shaping this correlation. However, such an analysis is beyond the scope of this study and is left for future work.

While Figure 11 illustrates the significance of individual parameters in influencing $\Delta f^j_{\mathrm{eff}}$ , it does not quantify the extent to which these external drivers may overlap or exhibit degeneracy in their correlation coefficients. To quantify potential levels of degeneracy between different external drivers, we proceed by computing the cross-correlation between these quantities.

4.2.2. Corr(External driver, External driver)

Figure 12 depicts the correlation coefficient matrix, illustrating the relationships between external parameters and $\Delta f^j_{\mathrm{eff}}$ , as well as among themselves. The left panel corresponds to $f_{\mathrm{MM}}=0.05$ , while the right panel illustrates the case with $f_{\mathrm{MM}}=0.2$ .

Figure 12. Depiction of the correlation coefficient matrix illustrating the relationships between external parameters and $\Delta f^j_{\mathrm{eff}}$ , as well as among themselves. The left panel corresponds to $f_{\mathrm{MM}}=0.05$ , while the right panel illustrates the case with $f_{\mathrm{MM}}=0.2$ .

Positive and negative cross-correlations are reported, although the statistical significance of negative correlations is notably smaller than that of positive cases. Hence, our focus lies solely on positive correlations. Notably, $\overline{\mathrm{MR}}$ exhibits strong correlations with $\overline{D}$ , $z_{\mathrm{merg}}$ , and $\log{(\overline{\frac{S2}{S1}})}$ . This is attributed to higher mass-ratios potentially exerting stronger impacts on average impact parameter values. Additionally, the positive correlation with $z_{\mathrm{merg}}$ suggests that mergers with higher mass-ratios may tend to occur at higher redshifts. Furthermore, a positive correlation with $\log{(\overline{\frac{S2}{S1}})}$ indicates that galaxy mergers with higher mass-ratios may also correspond to higher spin fractions, although alignment is not necessarily guaranteed.

$\overline{D}$ is correlated with $z_{\mathrm{merg}}$ , hinting that galaxy mergers at higher redshifts may involve larger impact parameters. While this observation is intriguing, further justification will be pursued in future work with a larger galaxy sample.

Additionally, the number of mergers shows a positive correlation with $\log{(S1 \cdot S2)_{\mathrm{M}}}$ , suggesting that a larger number of galaxy mergers statistically facilitates the occurrence of more aligned systems. Moreover, there are positive correlations between $\log{(\overline{\frac{S2}{S1}})}$ and $\log{(S1 \cdot S2)_{\mathrm{M}}}$ , indicating that, on average, systems with higher spin fractions may exhibit higher alignments.

These observations reveal overlapping effects between key drivers, ultimately enhancing their correlation coefficients with $\Delta f^j_{\mathrm{eff}}$ . Consequently, it is imperative to be cautious when listing the key external parameters, as their impacts may be lower due to the aforementioned degeneracy with other parameters.