Describing the Nonuniversal Galaxy Merger Timescales in IllustrisTNG: Effects of Host Halo Mass, Baryons, and Sample Selection

The IllustrisTNG simulations are a suite of magnetohydrodynamic cosmological simulations (F. Marinacci et al. 2018; J. P. Naiman et al. 2018; D. Nelson et al. 2018; A. Pillepich et al. 2018a; V. Springel et al. 2018; D. Nelson et al. 2019). These simulations are conducted using the moving-mesh code AREPO (V. Springel 2010) and incorporate a range of baryonic processes modeled as subgrid physics (R. Weinberger et al. 2017; A. Pillepich et al. 2018b). The cosmological parameters adopted in the simulations align with the Planck15 results (Planck Collaboration et al. 2016): Ω_m = 0.3089, Ω_Λ = 0.6911, Ω_b = 0.0486, h = 0.6774, σ₈ = 0.8159, and n_s = 0.9667. Among the suite of simulations, we utilize the TNG100-1 and TNG300-1 runs for this study. The TNG100-1 simulation features a box size of 75h⁻¹ Mpc, with 1820³ dark matter particles and 1820³ gas cells. This corresponds to mass resolutions of m_DM = 5.05 × 10⁶h⁻¹M_⊙ for dark matter and average m_gas = 9.44 × 10⁵h⁻¹M_⊙ for gas cells. The TNG300-1 simulation has a larger box size of 205h⁻¹ Mpc, containing 2500³ dark matter particles and 2500³ gas cells, with mass resolutions of m_DM = 3.98 × 10⁷h⁻¹M_⊙ for dark matter and average m_gas = 7.44 × 10⁶h⁻¹M_⊙ for gas cells.

Dark matter halos in the simulation were identified using the Friends-of-Friends method, while subhalos were cataloged using the SUBFIND algorithm (V. Springel et al. 2001). Subhalos across different snapshots were linked by merger trees constructed with the SUBLINK algorithm (V. Rodriguez-Gomez et al. 2015), which traces their formation history by assigning a unique descendant (if one exists) from the subsequent snapshot to each subhalo based on a scoring system that considers shared particles.

We use the DM-based SUBLINK merger tree to study the merger timescales of galaxies. We also tested the star- and gas-based merger trees and found minimal differences in the results. We define the infall time as the moment when the center of a subhalo crosses the virial radius r_vir of its host halo. The virial radius, r_vir, is the radius within which the mean mass density equals Δ_c (G. L. Bryan & M. L. Norman 1998) times the critical density of the Universe. In the SUBLINK merger tree, a merger is defined as the event in which two subhalos share the same descendant. For simplicity, we consider only mergers between a satellite subhalo and a central subhalo. Specifically, we include events where the descendant of the satellite subhalo lies in the main progenitor branch of the central subhalo at z = 0.

To ensure that both central and satellite galaxies are well resolved, we select host halos with virial masses M_vir(z = 0) > 10^10.7h⁻¹M_⊙ for TNG100-1 and M_vir(z = 0) > 10^11.7h⁻¹M_⊙ for TNG300-1 at z = 0, corresponding to more than 10,000 DM particles. Identifying satellite galaxies that have merged into central galaxies is nontrivial because of the challenges in robustly defining mergers in simulations. A merger, as defined by the SUBLINK merger tree, can correspond to one of three scenarios: a real merger, the actual disruption of a satellite (F. He et al. 2025), or the artificial disruption of a satellite (F. C. van den Bosch & G. Ogiya 2018). A real merger is typically defined as the point at which a large fraction of a satellite’s particles become mixed and indistinguishable from those of the central galaxy in phase space. In contrast, actual disruption refers to the physical disintegration of a satellite before its material merges with the central galaxy. The observed stellar streams in the Milky Way are remnants of such disruptions (K. Malhan et al. 2018). Artificial disruption, on the other hand, occurs when a satellite is entirely disrupted due to limited mass and force resolution; such satellites would remain intact in higher-resolution simulations. Distinguishing among these scenarios is essential for accurately studying galaxy merger timescales, since including artificial disruption events can significantly shorten the measured durations. Since distinguishing between real and artificial disruption is challenging, we first focus on genuine merger events. Results that include disruption cases will be discussed in Section 5.2. Therefore, following M. Boylan-Kolchin et al. (2008), we distinguish mergers from satellite disruptions using not only stellar mass or particle count but also the final specific angular momentum, j_merger, which, as shown by M. Boylan-Kolchin et al. (2008), is a better indicator than r_merger. In Figure 1, we show the joint distribution of the infall stellar mass of satellite galaxies and the ratio j_merger/j_infall of their specific angular momentum before the merger to that at infall. The stellar mass of galaxies is defined as the total stellar mass enclosed within twice the stellar half-mass radius. We show the results for both TNG100-1 and TNG300-1 for all systems with z_infall < 4. It is evident that j_merger/j_infall decreases with increasing stellar mass. Moreover, at a fixed stellar mass, TNG100-1 exhibits, on average, lower j_merger/j_infall values compared to TNG300-1. This suggests that galaxies with fewer particles are more susceptible to disruption, consistent with the results reported in F. C. van den Bosch & G. Ogiya (2018). To construct the merger sample, we apply both stellar mass and angular momentum cuts, as indicated by the red lines in Figure 1. Specifically, we use an infall stellar mass cut of for satellite galaxies in TNG100-1 and for those in TNG300-1, along with an angular momentum cut of j_merger/j_infall < 0.1. This stellar mass cut ensures that infalling subhalos contain more than 200 star particles and 5000 dark matter particles, guaranteeing that the subhalos are well resolved and that baryonic effects are present. After applying the selection criteria, we identify a total of 3792 mergers in TNG100-1 and 22,635 mergers in TNG300-1.

Zoom In Zoom Out Reset image size

Merger timescales refer to the time between a satellite’s infall into the host halo and its eventual disruption or merger with the central galaxy. The goal of this paper is to investigate the prediction of merger timescales for satellite galaxies based on their initial conditions at the time of infall. To achieve this, we calculate the infall and merger times, along with the orbital parameters and mass ratios at the infall time, for these satellite galaxies.

Since the simulations are stored in discrete snapshots, determining the exact infall and merger times requires approximations. For the infall time, we use the positions, velocities, and times of satellite galaxies from the snapshots immediately before and after infall. Assuming constant acceleration in the radial direction, we solve for the precise time at which the satellite galaxies first cross r_vir of their host halos. For the merger time, following C. Y. Jiang et al. (2008), we adopt the midpoint between the snapshots immediately before and after the merger as the merger time. The dense sampling provided by the large number of snapshots results in a time interval of 0.1–0.2 Gyr between snapshots, which only significantly affects the rare, very fast merger events.

For the initial conditions, all quantities are calculated at the snapshot immediately preceding infall. Three sets of parameters are typically required to determine the merger timescale: the mass ratio M_sat/M_host, the orbit parameters Θ, and the dynamical time T_dyn. For the mass ratio, we define M_host as the virial mass M_vir of the host halo and M_sat as the bound subhalo mass of the satellite. The dynamical time T_dyn is defined as , where r_host is the virial radius r_vir of the host halos. For a spherical potential, the orbit of a satellite can be fully determined by its energy and angular momentum (E, L). In practice, two dimensionless quantities are commonly used: and r_c(E)/r_host. Here,  = j_infall/j_c(E) represents the circularity of the orbit, while r_c(E) and j_c(E) denote the radius and specific angular momentum, respectively, of a circular orbit with the same energy E. As in many studies (C. Lacey & S. Cole 1993; C. Y. Jiang et al. 2008; L. Jiang et al. 2015), we adopt an SIS potential for halos to calculate the initial orbit parameters, as it is more convenient than the Navarro–Frenk–White (NFW; J. F. Navarro et al. 1997) potential. The differences in orbit parameters between the SIS and NFW potentials have been explored in L. Jiang et al. (2015, see their Figure 11). While the NFW potential may be more realistic, we find that the parameters derived from the SIS potential already provide accurate predictions for merger timescales, as indicated by below. Given its simplicity and effectiveness, we retain the use of the SIS potential in our analysis.

In Figure 2, we show the distributions of the initial conditions M_sat/M_host, , and r_c(E)/r_host for our merger samples. Despite the differences in the selection criteria for and M_vir between TNG100-1 and TNG300-1, their merger samples exhibit very similar initial distributions. Similar to C. Y. Jiang et al. (2008), most mergers in our sample have a mass ratio larger than 0.03, as smaller satellites cannot merge within a Hubble time. Additionally, we find that our merger sample exhibits smaller values compared to the complete infall sample studied in L. Jiang et al. (2015), while the r_c(E)/r_host distribution is quite similar. This arises because radial orbits have shorter merger timescales and are more likely to be included in our merger sample. The impact of these selection effects on our results will be examined in Section 5.2 using the full infall sample rather than only the merger sample.

Zoom In Zoom Out Reset image size

As in previous studies (C. Y. Jiang et al. 2008; V. Simha & S. Cole 2017; R. J. J. Poulton et al. 2021), we begin by attempting to construct a universal model for the entire merger sample. We adopt the exponential form of f() from M. Boylan-Kolchin et al. (2008), as it provides a better fit to our results than alternative forms such as power laws. The merger timescale model we adopt is

Here, α, β, and γ control the dependence of merger timescales on orbital circularity, mass ratio, and orbital energy, respectively, while A serves as a normalization factor.

We fit this model to our merger samples from TNG100-1 and TNG300-1 separately using the following expression for chi-squared:

where N is the number of mergers in the sample, p is the number of free parameters, represents the merger timescales calculated from simulations, and denotes the timescales predicted by the model. The χ² is minimized using the Levenberg–Marquardt algorithm. We obtain the following best-fit parameters for TNG100-1:

with a logarithmic scatter of . For TNG300-1, the best-fit parameters are

with a logarithmic scatter of .

In Figure 3, we compare the merger timescales from simulations and best-fit models. Each component in Equation (5) is examined separately, with the simulation data showing merger timescales normalized by the best-fit models that exclude the corresponding component. We find that the dependence of merger timescales on all components—, M_sat/M_host, and r_c/r_host—is well captured by Equation (5). The overall scatter in both TNG100-1 and TNG300-1 is small with smaller than 0.15, indicating that our model provides accurate predictions for merger timescales. However, although our model fits the merger timescales in both TNG100-1 and TNG300-1 samples well, the best-fit parameters differ, particularly for α, which describes the dependence on the orbit circularity . This suggests that galaxy merger timescales may not be universal and could depend on sample properties.

Zoom In Zoom Out Reset image size

The differences in best-fit models between TNG100-1 and TNG300-1, along with inconsistencies in previous studies such as C. Y. Jiang et al. (2008) and M. Boylan-Kolchin et al. (2008), suggest that galaxy merger timescales may depend on sample properties. Therefore, we divide the merger samples from TNG100-1 and TNG300-1 into different infall M_host bins and model them separately, as host halo mass is a fundamental property that influences key factors such as halo formation history, concentration, ellipticity, and subhalo distributions.

The best-fit parameters and their uncertainties are shown in Figure 4, while the data and best-fit models are presented in Figures A1 and A2, respectively. We find that the parameters indeed depend on M_host, and the results in TNG100-1 and TNG300-1 are now consistent. Specifically, the orbital parameters α and γ decrease with increasing M_host, while remaining relatively unchanged for M_host < 10¹²M_⊙. This suggests that the dependence of galaxy merger timescales on the initial orbit weakens in higher-mass host halos. For the mass ratio dependence, β increases slightly with M_host, from 0.75 to approximately 0.9. Additionally, the scatter also increases with M_host, further indicating a weaker connection between merger timescales and initial conditions in more massive host halos.

Zoom In Zoom Out Reset image size

Comparing our results with previous findings from C. Y. Jiang et al. (2008) and M. Boylan-Kolchin et al. (2008), we find that, for M_host < 10¹²M_⊙, the parameters α and γ approach the values of 1.9 and 1.0, respectively, as reported in the idealized simulations of M. Boylan-Kolchin et al. (2008). At higher host masses, our results align with the conclusion of C. Y. Jiang et al. (2008), which suggests that the dependence on initial orbital conditions becomes relatively weak. This may explain the discrepancies between these two studies as they are studying systems in different conditions. It is reasonable that in the idealized simulations of M. Boylan-Kolchin et al. (2008) merger timescales show the strongest connection to the initial conditions. In realistic simulations, violent processes in host halos, such as major mergers (N. E. Drakos et al. 2019a, 2019b) and interactions between satellite galaxies, can make the infall process more collisional, thereby erasing initial conditions and weakening their connection to merger timescales. For example, a sudden shift in the halo center position following a major merger can significantly weaken the connection, as it alters the reference frame used to calculate the orbital parameters. According to our results, these processes may occur more frequently in more massive host halos.

To characterize the dependence of merger timescales on M_host, we parameterize A, α, β, and γ as functions of M_host. Based on Figure 4, we adopt the following parameterization:

and a constant β = β₁, where is the error function, and is the complementary error function. These functional forms based on the error function are chosen to interpolate between two limiting cases: the strongest orbital dependence, as found in the idealized simulation by M. Boylan-Kolchin et al. (2008) with α = 1.9 and γ = 1.0, and the weakest case of no dependence, corresponding to α = 0 and γ = 0. Therefore, the model contains nine parameters: . We fit the model directly to the merger timescales of the combined TNG100-1 and TNG300-1 sample. The best-fit model is represented by black solid lines in Figure 4, with the corresponding parameters listed in the first row of Table 1. We find that the best-fit model generally follows the trends observed when modeling different M_host bins separately. Additionally, we obtain a similar as in the M_host bin results, indicating that our model provides a sufficiently accurate description.

Table 1. Descriptions and Best-fit Parameters of the Four Merger Timescale Models Constructed in This Study, Based on Different Samples, Starting Radii, and Model Parameterizations

Sample	Start	Model	Ah				αh		γh		β1	β2	Mβ
merger hydro	rvir	Equations (5), (7)	0.358	12.447	0.551	0.756	0.537	0.575	0.193	0.405	0.806	…	…
merger DMO	rvir	Equations (5), (7)	0.257	13.228	0.885	0.987	0.393	0.208	0.332	0.038	0.824	…	…
merger hydro	any	Equations (7), (8), (10), (11), (12)	0.360	12.376	0.605	0.833	0.406	0.700	-0.026	0.892	0.814	94.736	−0.456
infall hydro	any	Equations (7), (8), (10), (11), (12)	0.323	12.484	0.555	0.657	0.328	0.497	0.326	0.408	1.077	30.890	−0.026

Note. The first two rows use the merger sample starting at r_vir, from baryonic and DMO runs, respectively. The third uses the merger sample with arbitrary starting radii in baryonic runs. The fourth, recommended for general use, is based on the full infall sample with arbitrary starting radii in baryonic runs.

Download table as:  ASCIITypeset image

In Figure 4, we also present for different M_host using the parameters from the full-sample fit shown in Figure 3. The from the full-sample model is higher than that from the M_host-dependent model, confirming that incorporating M_host dependence improves the model’s accuracy. We also note that TNG300-1 shows a smaller than TNG100-1 at fixed M_host. This is because the merger samples in TNG100-1 span a broader mass ratio range—extending to lower values—due to a lower cut. If merger timescales retain some dependence on M_sat, this broader range would naturally lead to a larger in TNG100-1, which appears to be the case. We do not pursue a detailed investigation of the M_sat dependence, as the current model already yields a relatively low .

To investigate the differences between DMO and hydrodynamic simulations, we use the DMO simulations TNG100-1-Dark and TNG300-1-Dark, which serve as the DMO counterparts to TNG100-1 and TNG300-1. These simulations share the same cosmology, box sizes, and initial conditions. The mass resolutions of TNG100-1-Dark and TNG300-1-Dark are m_DM = 6.00 × 10⁶h⁻¹M_⊙ and 4.73 × 10⁷h⁻¹M_⊙, respectively. IllustrisTNG provides matched catalogs between DMO and baryonic runs using the SUBLINK weighting algorithm, which is similar to the method used for tracking descendants across different snapshots within the same simulation.

We start with a statistical comparison between DMO and baryonic runs, calculating merger timescales and their dependence on initial conditions in the DMO simulations using the same methodology as in the baryonic simulations. To construct the merger sample, we apply a selection similar to that in Figure 1. Specifically, we select host halos with virial masses of M_vir > 10^10.7h⁻¹M_⊙ for TNG100-1-Dark and M_vir > 10^11.7h⁻¹M_⊙ for TNG300-1-Dark at z = 0. Additionally, we require subhalo bound masses at the infall snapshot to be M_sat > 10^10.5h⁻¹M_⊙ for TNG100-1-Dark and M_sat > 10^11.5h⁻¹M_⊙ for TNG300-1-Dark. The same selection criteria of j_merger/j_infall < 0.1 and z_infall < 4 are also applied. Next, we calculate the merger timescales and initial conditions at the infall time and derive the merger timescale models for different infall host halo masses M_host according to Equation (5). Similar to the baryonic runs, we fit the results from the DMO merger sample using the models in Equation (7). The best-fit parameters are provided in the second row of Table 1, and the corresponding best-fit models are shown in Figure 5.

Zoom In Zoom Out Reset image size

Figure 5 compares merger timescale parameters between DMO and baryonic runs. We find that α is lower in DMO runs, with the difference increasing as M_host decreases, while β and γ remain largely similar. As a result of the variation in α, the normalization parameter A also differs. To investigate the cause of this difference, we compare the best-fit models across different M_host bins for TNG100-1-Dark and TNG100-1 in Figure A3, as well as for TNG300-1-Dark and TNG300-1 in Figure A4. We find that, in low-mass host halos, merger timescales in DMO simulations are generally longer than those in hydrodynamic simulations for highly radial orbits with smaller . To determine what is the cause of the difference, we further conduct a detailed one-to-one comparison of merger events in DMO and baryonic runs.

Since IllustrisTNG provides matched catalogs between DMO and baryonic runs, it offers a more effective way to study the factors influencing merger timescales through a one-to-one comparison. To match the events, we first identify the host halos in the DMO runs that correspond to those in the baryonic runs at z = 0. Then, for each merged satellite in the baryonic runs, we search for its DMO counterpart among the progenitors of the matched host halos at the satellite’s infall snapshot. Using this approach, we identify 3103 matched pairs out of 3792 in TNG100-1-Dark and 18,871 out of 22,635 in TNG300-1-Dark.

In the left panel of Figure 6, we compare the merger timescales T_baryon from the baryonic runs with T_dark from the corresponding matched events in the DMO runs. The results exhibit a quantized pattern because merger times are determined using the midpoint between the snapshots immediately before and after the merger. The differences between baryonic and DMO runs appear larger for short-duration mergers, while they remain fairly consistent for long-duration mergers. This is consistent with the results of the statistical comparison, as short-duration mergers typically occur in more radial orbits with smaller . To better quantify the dependence of the merger timescale difference on M_host and , the middle panel presents the logarithmic difference across different M_host and bins. For each event, we choose M_host and from the baryonic runs. We find that events with smaller indeed exhibit larger differences. However, interestingly, events with large and small M_host exhibit different trends. For M_host < 10^12.5h⁻¹M_⊙, merger timescales are larger in the DMO runs, whereas, for M_host > 10^12.5h⁻¹M_⊙, they are larger in the baryonic runs. Additionally, the magnitude of this difference increases toward both smaller and larger M_host.

Zoom In Zoom Out Reset image size

The differing impact of baryons on either side of M_host ∼ 10^12.5h⁻¹M_⊙ suggests that baryonic effects can influence merger timescales in multiple ways, either extending or shortening them. The dominant effect may vary depending on M_host, leading to the observed trend. First, the presence of a disk structure in the central galaxy can more effectively disrupt infalling satellites, shortening merger timescales—particularly for those on radial orbits that pass through the central region (G. Toth & J. P. Ostriker 1992; H. Velazquez & S. D. M. White 1999; A. Helmi et al. 2003; S. Garrison-Kimmel et al. 2017; J. Samuel et al. 2020). Second, the presence of baryons in the satellite galaxy may increase the inner density of the subhalo, making it more resistant to tidal disruption. If the tidal forces in the inner region of the host halo are strong enough to disrupt the subhalos in DMO simulations, the survival of the corresponding subhalo in baryonic simulations may lead to longer merger timescales. However, if the tidal forces are not strong enough to completely disrupt the subhalo but instead significantly reduce its mass in the DMO simulation, longer merger timescales in DMO runs are expected, as reported in M. Boylan-Kolchin et al. (2008). Finally, in halos hosting high-mass elliptical galaxies, strong feedback from the supermassive black hole reduces the matter density in the central region of the host halos (K. Xu et al. 2024, see their Figure 17), primarily within r < 100h⁻¹ kpc. This reduction in density may weaken the deceleration caused by dynamical friction, as described by Chandrasekhar’s formula, leading to an increase in merger timescales.

Since one of the most significant differences between low- and high-mass galaxies is their morphology—where high-mass galaxies are mostly elliptical while low-mass galaxies are predominantly disk shaped—we suspect that the difference in merger timescales for low-mass halos is due to the impact of their disk structure. To support this hypothesis, the right panel of Figure 6 presents the logarithmic difference in merger timescales, , across different M_host and F_disk bins, where F_disk represents the disk-to-total stellar mass fraction of the central galaxy at the merger snapshot (S. Genel et al. 2015). Disk particles are identified based on the criterion J_z/J(E) > 0.7, where J_z is the specific angular momentum of the particle, and J(E) is the expected angular momentum for a circular orbit at the same location. We find that depends on F_disk, with larger values observed for higher F_disk, partially supporting our hypothesis. For high-mass halos, we believe that the resistance of baryons to tidal disruption may be the dominant effect, as many subhalos are disrupted much more quickly in DMO runs. This is evident in the left panel of Figure 6, where numerous events exhibit in the DMO runs but in the baryonic simulations.

We conduct a more detailed analysis by comparing the infall history of matched merger events in DMO and baryonic runs. Figure 7 presents two matched events, one in a low-mass host halo hosting a central disk galaxy and the other in a high-mass host halo with a central elliptical galaxy. These two events have similar initial conditions, with  ∼ 0.13 and M_sat/M_cen ∼ 0.1. The selected events have relatively radial orbits, which exhibit more pronounced differences between the two runs. The mass ratios are chosen to ensure sufficiently long merger histories for clearer illustration, while avoiding rare cases with very small mass ratios. Notably, the differences between the DMO and baryonic runs exhibit opposite trends in high-mass and low-mass halos. In the low-mass halo, the subhalo in the baryonic run loses a significant amount of mass after its first pericenter passage and merges with the central galaxy immediately after the second passage. In contrast, in the DMO run, the subhalo survives through four pericenter passages. In the high-mass halo, the subhalo in the DMO run experiences a more rapid loss of mass and sinks to the center faster than in the baryonic run.

Zoom In Zoom Out Reset image size

While this analysis supports our hypothesis, further controlled experiments are necessary for a more comprehensive understanding, given the highly chaotic nature of the processes and the numerous influencing factors. Additionally, a thorough investigation is needed to determine whether these effects are physical or merely numerical artifacts.

Since dark matter halos exhibit universal density profiles (J. F. Navarro et al. 1997), a natural question arises: can our model be extended to merger events initiated at any radius within the halo, as the inner regions of halos can be considered as new halos with different overdensities (accounted by T_dy) and concentrations? Some studies have claimed a validation for their models (V. Simha & S. Cole 2017; R. J. J. Poulton et al. 2021). Therefore, we extend our model in Equation (5) to estimate merger timescales for any initial starting radius r as follows:

where M_host(<r) is the mass of the halo enclosed within the radius r, is the dynamical time, M_sat represents the bound mass of the subhalo, and the orbital parameters and r_c are still computed using the SIS potential for convenience, constrained by (r, M_host(<r)). We first maintain the mass ratio dependence as g(x) = x^β. The mass dependence of α, β, and γ is preserved as a function of the virial mass at the first crossing of r_vir, following Equation (7), with the best-fit values given in the first row of Table 1.

To test the model, we use the merger sample from TNG100-1 and TNG300-1. For each event, we treat each intermediate snapshot during its infall history as the starting point of a new event and compute the corresponding initial conditions and merger timescales. This process requires particle data from each snapshot to determine M_host(<r). In Figure 8, we present the logarithmic difference, , between simulation results and model predictions from Equation (8). This difference is shown as a function of , mass ratio M_sat/M_host(<r), orbit energy r_c/r, host halo mass M_host,0 at the time of crossing r_host,0 = r_vir, and starting radius r. We find that, after entering the halo (), our extended model underestimates the merger timescales, suggesting a transition before and after infall. Moreover, our model significantly overestimates the merger timescales for mass ratios <10⁻², where fewer cases were present in the previous sample used for parameter constraints. This may be reasonable, as subhalos in low mass ratio events experience weaker dynamical friction but are more susceptible to disruption. As a result, this effect may cause events with low mass ratio to deviate from the previous mass dependence.

Zoom In Zoom Out Reset image size

To better describe the merger timescales for events starting at any radius within the halo, we introduce a correction based on r and adjust the mass ratio dependence at the low-mass end. We find that the dependence of on r originates from M_host(<r). Figure 9 shows M_host(<r)/M_host,0 as a function of r/r_host,0 for our merger sample. The fraction M_host(<r)/M_host,0 increases with r/r_host,0, following a power-law trend up to r/r_host,0 ∼ 10^−0.5, beyond which the growth gradually slows down. This is due to the halo’s nonspherical mass distribution and boundary. As r approaches r_vir, the enclosed sphere extends beyond the irregular halo boundary, capturing little additional mass and slowing the total mass increase. If M_host(<r) determines the merger timescale, this is not an issue. However, if the relevant factor is not the enclosed spherical mass but an alternative effective mass, such as the elliptical mass, this trend could introduce a bias across different r, as demonstrated in Figure 8. This bias can be corrected by defining an effective host halo mass. The relationship between M_host(<r)/M_host,0 and r/r_host,0 is well described by

with best-fit parameters r_t = 0.833 and η = 1.257. Based on this relation, we define a correction factor

which is normalized such that f_c(1) = 1, ensuring consistency with our previous results. The effective host halo mass is then defined as f_cM_host(<r), which is also presented in Figure 9. To better capture the mass ratio dependence, motivated by the power-law upturn bias seen in Figure 8, we replace with a double power-law form:

where M_ratio = M_host(<r)/M_sat. We test replacing M_host(<r) with f_cM_host(<r) in M_ratio, Λ, and T_dy and find that only T_dy requires modification. Adjusting the other components has negligible or adverse effects on the results. Therefore, we define a new T_dy as

Zoom In Zoom Out Reset image size

By combining Equations (7), (8), (10), (11), and (12), we develop a model to describe merger timescales for mergers starting at arbitrary starting radii. This model includes 11 parameters: . When including r_t and η, the total number of parameters increases to 13. We fit this new model to the merger sample starting at arbitrary starting radii. The best-fit parameters are listed in the third row of Table 1, and the residuals as a function of various parameters are shown in Figure 10. Our new model effectively describes the merger timescales for the arbitrary starting radius sample, achieving , with the scatter showing little to no dependence on any parameters. In Figure 10, we also show the scatter for the best-fit model without the correction for T_dy in the r/r_host,0 panel. Including this correction really helps reduce the bias dependence on r/r_host,0. In the new model, we find a weaker dependence of α and γ on M_host,0 compared to the previous model, with α_h = 0.406 (previously 0.537) and γ_h = −0.026 (previously 0.193). Notably, γ now shows almost no dependence. This suggests that M_host,0 may not be a reliable indicator of system variation for mergers that do not start at r_vir, and more fundamental parameters might be needed. However, since our model already performs well and M_host,0 remains the most accessible quantity, we do not explore alternative options.

Zoom In Zoom Out Reset image size

One remaining concern regarding our merger timescale results is selection effects, a common issue in studies of merger timescales using cosmological simulations. Since we only analyze satellites that have already merged with centrals, our sample is likely biased toward events with shorter merger timescales. If this sample is taken as a representative of all infall events, our model may underestimate the true average merger timescales. With the model for arbitrary starting radii, we can now use the entire infall sample—whether the satellites have merged, have been disrupted, or are still surviving—to study selection effects. This is possible because our model allows us to predict merger timescales for any two points in an event’s timeline and compare the predicted time intervals with those from simulations. As a result, satellites do not need to have already merged with the central to be included in the analysis.

To construct the infall sample, we select all satellites that have crossed r_vir after z = 4, using only the mass selection criterion as before: M_vir(z = 0) > 10^10.7h⁻¹M_⊙ and for TNG100-1, and M_vir(z = 0) > 10^11.7h⁻¹M_⊙ and for TNG300-1. In Figure 11, we present the distributions of the initial properties—, , and r_c/r_host—at the time of first crossing r_vir for the infall sample. We also show the contributions from three subsamples: satellites that have merged with the centrals (), those that have been disrupted (, regardless of whether the disruption is physical or artificial), and those that are still surviving at z = 0. If we use only the merger sample instead of the full infall sample as before, we tend to select satellites with higher mass ratios, more radial orbits, and possibly shorter merger timescale events for the same initial conditions, which may bias our model toward underestimating merger timescales.

Zoom In Zoom Out Reset image size

Therefore, we refine our model by incorporating the entire infall sample. For the merger sample, we maintain the previous approach, requiring only the prediction of the merger timescale at the infall time, T_merger(t_infall), where t_infall is defined as the time of first crossing r_vir. For the disruption and surviving samples, we predict both T_merger(t_infall) and T_merger(t₅₀), where t₅₀ is the last time the subhalo contained more than 50 particles. If the subhalo still has more than 50 particles at z = 0, then t₅₀ corresponds to the present time. The requirement of at least 50 remaining particles is imposed to reduce numerical artifacts. With these predictions, we compute ΔT_merger = T_merger(t_infall) − T_merger(t₅₀) from our model and compare it with simulation results. For the merger sample, ΔT_merger simply equals T_merger(t_infall). We have also tested t₁₀₀ and found similar results, confirming the robustness of this choice.

We first assess whether our previously derived arbitrary starting radii model, based on the merger sample, accurately predicts ΔT_merger for the entire infall sample. Figure 12 compares model predictions of ΔT_merger with simulation results. Overall, the model performs well across all three subsamples, except when the mass ratio is very small (M_sat/M_host < 10⁻³). For the merger sample, T_merger from simulations is generally slightly shorter than the model predictions, with discrepancies reaching 0.2 dex for M_sat/M_host < 10⁻². This aligns with expectations, as the merger sample tends to favor shorter merger timescales. Despite being calibrated on the merger sample, the model’s ability to incorporate information from arbitrary radii appears to mitigate selection effects, likely due to the collisional nature of the merger process.

Zoom In Zoom Out Reset image size

We refine our model by fitting ΔT_merger from the entire infall sample, using the same framework as the arbitrary starting radii model. The best-fit parameters are listed in the fourth row of Table 1, and Figure 13 compares model predictions with simulation results. The new model shows improved accuracy, particularly for low mass ratios (M_sat/M_host < 10⁻³) and high r_c/r, with an average logarithmic scatter of for the entire infall sample. For M_sat/M_host < 10⁻², T_merger from the merger sample remains 0.2 dex lower than the model prediction, highlighting the presence of selection effects. This also indicates that incorporating the entire infall sample into the fitting framework effectively reduces selection effects. However, this comes at the cost of increased scatter due to the need for predictions at two different times.

Zoom In Zoom Out Reset image size