Calibrating Galaxy Infall Times in Groups and Clusters with IllustrisTNG Simulations

IllustrisTNG is a suite of cosmological, magnetohydrodynamical simulations of galaxy formation (D. Nelson et al. 2017, 2019; A. Pillepich et al. 2017; V. Springel et al. 2017; J. P. Naiman et al. 2018; F. Marinacci et al. 2018). In this work we use two runs from this suite, TNG300 and TNG-Cluster:

1. 1.  
   TNG300. A large-volume simulation that evolves a (302.6 Mpc)³ cube from z = 127 to z = 0.
2. 2.  
   TNG-Cluster. An extension of the IllustrisTNG project containing 352 high-resolution zoom-in re-simulations of galaxy clusters (D. Nelson et al. 2024), evolved from z = 120 to z = 0.

Both simulations adopt the same physical model, and the zoom-in targets in TNG-Cluster have the same numerical resolution as TNG300 (m_baryon ≈ 1.1 × 10⁷ M_⊙, m_DM ≈ 5.9 × 10⁷ M_⊙). The regions between zoom-in targets in TNG-Cluster have lower resolution.

TNG300 contains many group- and cluster-mass halos but relatively few very massive clusters (M₂₀₀ > 10^14.5 M_⊙). Combining TNG300 with TNG-Cluster therefore yields a large, statistically powerful sample spanning a wide halo mass range (M₂₀₀ ∼ 10¹³–10^15.5 M_⊙).

The simulations provide 100 snapshots with associated group catalogs, identifying halos (groups/clusters) via a friends-of-friends algorithm and subhalos (galaxies) via the SUBFIND algorithm (V. Springel et al. 2001). Merger trees are constructed using SubLink (V. Rodriguez-Gomez et al. 2015) and LHaloTree (V. Springel et al. 2005).

2.1.1. Sample

For each subhalo, we extract the position, velocity, and stellar mass. We remove subhalos flagged as no longer physically meaningful (e.g., tidally disrupted or stripped below the resolution limit; D. Nelson et al. 2019). For each halo, we extract the central subhalo ID and virial mass (M_h), assigning the position and velocity of the central galaxy to those of its host halo. In TNG-Cluster, we keep only the halos corresponding to the 352 primary zoom-in targets, to ensure consistent resolution with TNG300. In TNG300, we keep all 3545 groups (10¹³ M_⊙ < M_h < 10¹⁴ M_⊙) and 280 clusters (10¹⁴ M_⊙ < M_h < 10^15.2 M_⊙) from the catalog.

Our galaxy sample includes all subhalos with stellar masses >10⁹ M_⊙, ensuring each galaxy contains at least 90 particles (sufficient for this work, where we are interested in where the subhalos are in the halo, not their detailed resolved properties) and is well-matched to observational surveys (V. M. Sampaio et al. 2021; M. Oxland et al. 2024). Throughout this paper, we refer to subhalos as “galaxies.” We restrict our analysis to galaxies associated with groups and clusters of mass M_h > 10¹³ M_⊙, and that lie within three virial radii (3R₂₀₀) of the host center. This selection is based on a geometric criterion rather than on friends-of-friends membership, allowing the calibration to be more directly applicable to observational samples constructed with different group-finding algorithms.

In order to mimic observations, we select a line of sight (LOS) and project position data along this LOS, retaining only two spatial coordinates and one velocity coordinate (v_LOS). We use the three principal axes (x, y, z) as independent LOS projections for both simulations to increase statistics.

Our final sample contains more than 1.5 × 10⁶ galaxies (about 521,000 per line of sight), of which roughly 520,000 are members of groups or clusters. The other galaxies are either field galaxies or members of groups with M_h < 10¹³ M_⊙. The distributions of host halo masses and galaxy stellar masses in our final sample are presented in Figure 1. The top panel shows that combining the two simulations greatly increases the halo mass range we have access to.

Zoom In Zoom Out Reset image size

The definition of infall time varies in the literature. Here, we define a galaxy as having entered a cluster or group, and thus become a satellite of this halo, when it first crosses three times the virial radius (R₂₀₀) of its current host. The corresponding infall time () is the time elapsed between this crossing and the present day, for the case of z = 0.

While some studies define infall time as the moment a galaxy first crosses R₂₀₀, we instead adopt 3R₂₀₀ for two reasons. First, environmental effects (e.g., ram pressure stripping, tidal forces) are known to influence galaxies beyond R₂₀₀ (M. L. Balogh et al. 2000; C. P. Haines et al. 2015; K. A. Oman et al. 2020). Second, because our analysis examines PPS out to 3R₂₀₀, adopting R₂₀₀ would introduce an artificial discontinuity in the distribution at R₂₀₀ in PPS, complicating model fitting. This choice is discussed further in Section 3.3.

We measure using SubLink merger trees from the IllustrisTNG database, which provide galaxy positions and the time evolution of host halo radii. For each galaxy, we compute its 3D distance from the central galaxy. If the galaxy lies within 3R₂₀₀ of its host at the reference redshift (e.g., z = 0), we then trace its normalized group-centric distance, R_3D(z)/R₂₀₀(z), backward through all earlier snapshots. The earliest snapshot at which this ratio first falls below 3 defines the infall redshift, and the corresponding lookback time from the reference redshift gives the galaxy’s .

We construct PPS diagrams including all satellite galaxies, where PPS coordinates are defined as

where σ is the halo velocity dispersion, computed following X. Yang et al. (2007):

Figure 2 compares the distributions obtained when defining infall at 3R₂₀₀ (top) versus R₂₀₀ (bottom). Using R₂₀₀ introduces a clear discontinuity, at r = 1 in PPS, making interpretation difficult. The 3R₂₀₀ definition avoids this issue and includes galaxies that are already environmentally affected at large radii. In both panels, galaxies are color-coded by the mean within each PPS bin.

Zoom In Zoom Out Reset image size

Due to projection effects, the PPS region shown in Figure 2 is contaminated by galaxies located beyond 3R₂₀₀ from the halo center in three-dimensional space. We classify as interlopers those galaxies whose true 3D distance from the halo center lies between 3R₂₀₀ and 10R₂₀₀ (following H. Dou & H. Yu 2025) but whose projected position falls within 3R₂₀₀. Such galaxies can mimic genuine satellites in PPS despite never having experienced comparable environmental conditions, potentially biasing inferred infall times and quenching trends. All interlopers were removed prior to constructing the distributions shown in Figure 2.

Figure 3 shows how the interloper fraction varies across PPS. As found in previous studies (e.g., J. Rhee et al. 2017; H. Dou & H. Yu 2025), contamination is highest at large radii (r ≳ 2) and high relative velocities, where interlopers can account for more than 60% of apparent satellites. In the central regions (r ≲ 1), the fraction falls below 10%. We exclude interlopers (i.e., galaxies with true 3D distance >3R₂₀₀ but projected within 3R₂₀₀) from our analysis, since we have no reliable estimate of their infall times. Observationally, contamination increases sharply at r > 2; we therefore recommend restricting analyses to galaxies within r < 2.

Zoom In Zoom Out Reset image size

Infall time varies systematically with position in PPS, but each location in PPS exhibits substantial intrinsic scatter in . For practical use with spectroscopic samples, we seek a continuous mapping that captures the dominant dependences with as few parameters as possible. We begin at z = 0 to identify an analytic relation between and a small set of predictors accessible to observers. Specifically, we consider the predictors , where M_∗ is stellar mass and M_h is host halo mass. Although previous studies (e.g., J. Rhee et al. 2017; A. Pasquali et al. 2019) showed that the dependences of infall time on stellar and host halo masses are weak, we included these parameters to test whether even small trends could improve the predictive performance of the calibration.

We model

and determine f in two stages. First, we coarsely grid the four-dimensional space and use symbolic regression to identify a compact functional form; second, we refit the coefficients of the best-fitting function using all individual galaxies (with interlopers removed).

We use the PySR symbolic regression package (M. Cranmer 2023) to search for functional forms relating to the four parameters. Symbolic regression provides an interpretable, data-driven way to uncover analytic relationships between projected phase-space coordinates and infall time, allowing us to capture the principal trends without imposing a predefined model. To keep the computation tractable, we sample the 4D space on regular grids, replacing all galaxies in each cell with a single point defined by the mean parameter values and mean infall time of that cell. Grid sizes range from 20 × 20 × 8 × 8 down to 8 × 8 × 6 × 1 cells (dimensions listed in the order r, v, , ). In this analysis, each cell is weighted by , where N_gal is the number of galaxies in the cell, so that sparsely populated regions have less impact on the fit. Results are consistent without weighting, albeit with slightly larger scatter.

To control model complexity and avoid overfitting, we restricted the symbolic-regression search space to the monotonic operators and prohibited highly nonlinear expressions (e.g., r^v). We repeated the search across multiple grid resolutions. Across these trials, only a small minority of candidate functional forms retained a dependence on halo mass. As described in the next section, comparison among all candidate models showed that those including M_h do not provide a better fit to the data than M_h-independent alternatives. A complementary correlation analysis of the full, unbinned sample likewise indicated that M_h adds negligible explanatory power once (r, v, M_∗) are included. Consequently, we removed M_h from the calibration: subsequent grid searches were performed without a halo mass dimension, and all remaining fits and results in this paper are based on (r, v, M_∗) only.

4.1.1. Best-fit Functional Form

We rank candidate expressions by their rms error (RMSE, ) and select the lowest-RMSE form at z = 0 returned by PySR:

We clip nonphysical negative predictions to zero:

After identifying the functional form using the gridded data, the coefficients were fitted using the full (interloper-free) galaxy sample, yielding

We note that this functional form has trends with r and v consistent with qualitative expectations: infall time increases toward smaller projected radii and smaller velocities, corresponding to galaxies that are deeper within their host potential. These trends are recovered in Equation (3). We emphasize that this functional form should be interpreted as a calibrated estimator rather than a physical model.

4.1.2. Performance at z = 0

Figure 4 compares the predicted from Equation (3) with the true infall times measured directly from the simulation at z = 0. In the top panel (purple points), each point shows the mean within a grid cell (grid resolution 20 × 20 × 6 in (r, v, M_∗)), while the predicted values are computed using coefficients refit on individual galaxies. The bottom panel presents the same comparison for individual galaxies. In both cases, the model captures the overall trend but compresses the dynamic range: short are biased high and long are biased low. This compression is particularly evident in the two-dimensional histogram of individual galaxies (bottom).

Zoom In Zoom Out Reset image size

For individual-galaxy predictions (bottom panel), with the A, B, and C values given above, the overall RMSE is 2.31 Gyr across the PPS. As shown in Figure 5, absolute errors are smallest in the upper right region of PPS, where typical values approach zero; however, the fractional error can be large (RMSE, see Appendix A). Conversely, absolute errors increase toward the lower left region, but they represent a smaller fraction of the true value because can reach up to 13.3 Gyr.

Zoom In Zoom Out Reset image size

At higher redshifts, less cosmic time has elapsed and clusters are still assembling, so the infall-time calibration may evolve with redshift. We first test whether the same functional form (Equation (3)) remains valid at earlier epochs by examining snapshots at z ∈ {0, 0.2, 0.3, 0.5, 0.7, 1}. To test this, we refit the coefficients A(z), B(z), and C(z) in Equation (3) at each redshift using curve_fit, and compute alternative candidate functional forms at each redshift following the method described in Section 4.1. We then compare the RMSE of infall times predicted by Equation (3) with the refitted coefficients to those of the alternative fits. We find that Equation (3) continues to provide a good description up to z = 1, with RMSEs similar to or lower than the alternatives.

Quadratic functions of z describe the evolution of A and B, while a linear function suffices for C (Figure 6). Incorporating these redshift-dependent coefficients reduces the RMSE compared to using a fixed z = 0 calibration at all redshifts.

Zoom In Zoom Out Reset image size

The redshift dependence is

Figure 6 shows the fitted coefficients for TNG300 and TNG-Cluster separately, and for the combined dataset. Although the absolute values differ slightly between simulations, the redshift trends are consistent, and the RMSE computed using the combined or single-simulation coefficients differs by less than 0.06 Gyr at each redshift. The increasing difference between the coefficients fitted to TNG300 and TNG-Cluster at high redshift, especially for coefficient B (Figure 6), could indicate a modest increase in sensitivity to halo mass at earlier epochs. However, because these differences have a negligible impact on the RMSE, we do not pursue this possibility further here.

Figure 7 shows the evolution of the mean RMSE across the PPS with redshift. The RMSE decreases steadily with increasing redshift, reflecting the shorter cosmic times available for infall at earlier epochs. Averaged over the entire PPS region, uncertainties are 2.5 Gyr at z = 0 and decline to 1 Gyr at z = 1. The RMSE evolution with redshift is well described by

Zoom In Zoom Out Reset image size

Variations in RMSE across PPS at each redshift are shown in Appendix A, Figure 11.

If the z = 0 coefficients are applied at higher redshift without refitting, the RMSE systematically increases with redshift, as shown by the orange triangles in Figure 7. These results underscore the importance of incorporating explicit redshift dependence into the calibration.

A. Pasquali et al. (2019, hereafter P19) divided projected phase space into eight zones out to R₂₀₀ and calibrated mean infall times within each region using simulations. Their definition of infall time corresponds to a galaxy’s first crossing of R₂₀₀, whereas our calibration adopts 3R₂₀₀ as the entry radius. Because of this difference, the absolute values are not directly comparable, but the relative trends can be examined.

Table 1 summarizes the mean predicted by our model within the P19 zones (excluding galaxies with r > 1 for consistency). The progression of mean across zones closely matches that reported by P19. The mean difference of 2.3–3.3 Gyr between the two definitions corresponds roughly to the time required for a galaxy to travel from 3R₂₀₀ to R₂₀₀ on first infall.

Table 1. Mean Infall Times in the P19 Zones and RMSE of Our Calibration as a Function of Redshift

Zonea	b	c	ΔTd	(Gyr)e
	(Gyr)	(Gyr)	(Gyr)	P19	z = 0.0	z = 0.2	z = 0.3	z = 0.5	z = 0.7	z = 1.0
1	8.41	5.42	2.99	2.51	2.30	1.82	1.66	1.35	1.12	0.95
2	7.67	5.18	2.49	2.60	2.43	1.93	1.76	1.46	1.23	1.01
3	6.99	4.50	2.49	2.57	2.44	1.92	1.74	1.44	1.24	1.04
4	6.20	3.89	2.31	2.34	2.37	1.83	1.63	1.33	1.14	0.95
5	5.70	3.36	2.34	2.36	2.45	1.87	1.67	1.35	1.13	0.88
6	5.38	2.77	2.61	2.29	2.47	1.90	1.68	1.33	1.13	0.88
7	5.12	2.24	2.88	1.97	2.40	1.82	1.61	1.29	1.07	0.84
8	4.69	1.42	3.27	1.49	2.12	1.60	1.40	1.18	1.09	0.95

Notes. ^aZone number from P19. ^bMean infall time computed with Equation (3) at z = 0. ^cMean infall time from P19. ^dDifference between our mean infall time and that of P19. ^eRMSE of predicted infall times from P19 and from our calibration at each redshift.

Download table as:  ASCIITypeset image

We compare the RMSEs in between our calibration and that of P19 at z = 0 by evaluating Equation (8) within the zones defined by P19. The z = 0 RMSEs from our model and those reported by P19 (fifth column) are listed in Table 1. Our RMSEs are similar to P19’s in the inner zones and somewhat larger in zones 7 and 8. This difference may reflect the weighting adopted in our calibration, which places greater statistical emphasis on the more densely populated inner regions of PPS. Alternative weighting schemes were explored, but none produced a lower global RMSE. Because our calibration defines infall at 3R₂₀₀ rather than R₂₀₀, the corresponding values are systematically higher with our method. Although our fit includes stellar mass, the overall RMSE values remain comparable to those of P19. The remaining scatter highlights the intrinsic difficulty of reducing the complex infall-time distributions to a single representative value per zone. For completeness, Table 1 also reports our RMSEs at higher redshifts, which are not available in P19.

Although our calibration does not yield lower RMSE than P19 at z = 0, its continuous functional form eliminates edge effects associated with discrete zone boundaries. Our framework extends naturally to higher redshift, making it directly applicable to other observational samples.

As in Section 4.1, we divide the three-dimensional space into regular cells, using the same range of grid sizes as before. In each cell, we compute the distribution of infall times and fit it as the sum of two Gaussian components using the GaussianMixture implementation from the scikit-learn library (F. Pedregosa et al. 2011). An example for a 6 × 6 × 1 grid is shown in Figure 8. The values of the Gaussian parameters used to plot the black curves of Figure 8 are listed in Appendix C, Table 3.

Zoom In Zoom Out Reset image size

In most PPS regions, the distributions are well described by a sum of two components. In cells with v ≳ 1.5, a single Gaussian component is often sufficient. Cells at very large radii and velocities (r + v > 4) contain relatively few galaxies and should be interpreted with caution.

We use PySR (M. Cranmer 2023) to derive symbolic expressions for the two Gaussian means, T_mean,1 and T_mean,2, as functions of (r, v, M_*). Note that T_mean,1 is always the component with the smaller . The best-fit expressions are

The average RMSE values across PPS, computed using a 20 × 20 × 6 grid, are RMSE₁ = 0.56 Gyr and RMSE₂ = 1.20 Gyr. We use gridded cells (20 × 20 × 6) here to allow measurement of full distributions within each region. Figure 9 shows that the measured Gaussian means and those predicted by Equations (10) and (11) agree well.

Zoom In Zoom Out Reset image size

The functional forms of Equations (10) and (11) differ substantially from that of Equation (8). They likely reflect a combination of projection effects in PPS and variations in accretion histories (e.g., early- versus late-infall populations), but we do not attempt to disentangle these contributions here. We emphasize that these relations should not be interpreted as a physical model of accretion. They provide an empirical decomposition of the infall-time distribution, illustrating that a single-valued estimator (Equation (8)) cannot fully represent what is intrinsically a multivalued distribution.

In principle, the full infall-time distribution at each PPS location could be predicted by fitting functional forms not only for the means but also for the standard deviations and relative weights of the two Gaussian components. Although the distributions appear bimodal and are well described by double Gaussians, our attempts to model the variation of these additional parameters across PPS were unsuccessful: all resulting models provided poor fits to the data, with RMSEs comparable to or larger than the intrinsic scatter of the data. This suggests that the variables considered here, (r, v, M_*), are insufficient to fully capture the shape of the underlying distributions.

Figure 10 shows the weight of the Gaussian with the smaller mean T_inf in PPS. A weight near unity indicates that the distribution is dominated by the first component. Note that regions with r > 2 or r + v > 4 are sparsely populated and more affected by interlopers in observational studies.

Zoom In Zoom Out Reset image size

A two-Gaussian description in PPS successfully captures distinct early- and late-infall populations and reduces the RMSE by roughly a factor of 2, provided the correct component is identified. However, the variables used here do not reliably predict which component an individual galaxy is in. In short, while the single-value calibration offers a practical, redshift-dependent tool for estimating infall times, explicitly incorporating the bimodality has the potential to further reduce scatter in future work. The next step is to identify additional physical parameters that can distinguish between the two populations and enable more precise reconstruction of environmental histories.