© The Authors 2026

1. Introduction

Galaxies form and evolve in dark matter haloes; therefore, we expect a strong correlation in the evolution of these two objects (White & Rees 1978). As haloes grow hierarchically, halo mergers become a fundamental way of increasing mass. While galaxies also merge, they do so on much longer timescales, so more massive haloes will contain several galaxies.

In a massive halo, the most massive and brightest galaxy is expected to be at the centre of the potential, with a velocity similar (but not identical) to that of the whole halo. This galaxy is commonly referred to as the central galaxy, whereas the remaining galaxies in the halo are known as satellite galaxies. When looking at the dark matter component of the halo, the position and velocity of satellites match those of subhaloes, which are the halo remnants of minor haloes that have fallen into the main halo. Using the Subhalo Abundance Matching formalism (SHAM, e.g. Vale & Ostriker 2006; Conroy et al. 2006), subhaloes can be associated with galaxies based on their positions, velocities, and rotational-velocity history. This technique, often used in the creation of mock galaxy catalogues, is capable of reproducing several galaxy clustering statistics, such as the two-point correlation function (in real and redshift space; Chaves-Montero et al. 2016; Ortega-Martinez et al. 2024), galaxy–galaxy lensing (e.g. Contreras et al. 2023a), and even higher-order statistics (e.g. Contreras et al. 2024), and has recently been used for the study of large-scale-structure effects such as assembly bias (e.g. Contreras et al. 2021a, 2023c), for characterising the clustering of the latest generation of galaxy surveys (e.g. Rocher et al. 2023; Ortega-Martinez et al. 2025), and even for deriving cosmological constraints (e.g. Contreras et al. 2023b; Mahony et al. 2026). Mocks using the SHAM formalism are normally run based on dark-matter-only simulations. The low computational cost of these simulations allows running them on larger volumes and in large numbers, which is necessary when analysing the observations provided by surveys such as DESI or EUCLID (e.g. Maksimova et al. 2021; Euclid Collaboration: Castander et al. 2025).

While the SHAM approach is accurate and computationally efficient, it has some known limitations. The presence of baryons, for instance, can modify the halo density profiles as well as the positions and velocities of subhaloes inside haloes. Differences between hydrodynamical and dark-matter-only orbital evolution have also been reported in merger studies (e.g. Pu et al. 2025). Because empirical models (such as SHAMs or halo occupation models) are usually calibrated on dark-matter-only simulations, and the positions and velocities of galaxies are set to be equal to those of subhaloes, dark matter particles, or to follow a standard NFW profile, the absence of baryonic effects will impact the clustering signal of these galaxies.

Modern abundance-matching models use additional criteria to determine which satellite subhaloes should be populated by a specific type of galaxy (Contreras et al. 2021b). These extensions are especially useful for accounting for selection effects, which refer to the fact that galaxy samples are built using different observational criteria. For instance, blue satellite galaxies are situated in the outer regions of haloes (Orsi & Angulo 2018), whereas redder galaxies are more prevalent in the inner regions. By adding this extra flexibility, models can populate subhaloes that better resemble the distribution of galaxies selected by a specific property (such as a cut in stellar mass or star formation rate). These mechanisms can also partially mimic the effects of baryons in the satellite distribution. While they do not shift the position of the dark matter subhalo as baryons do, they will populate subhaloes that are closer to the inner or outer part of the halo, producing a similar effect. However, to our knowledge, no comprehensive study has jointly quantified the impact of baryons on the positions and velocities of satellite galaxies and compared it directly to the impact of selection effects.

In this paper, we make use of a state-of-the-art cosmological hydrodynamic simulation, the MillenniumTNG simulation (MTNG, Pakmor et al. 2023; Hernández-Aguayo et al. 2023), to measure the impact of baryons on the positions and velocities of satellite subhaloes. By comparing the galaxy population of the hydrodynamic simulation (which models baryons and dark matter simultaneously) to that of a dark-matter-only simulation, run with the same initial conditions and cosmology as MTNG, we can determine, satellite by satellite, the shift in position and velocity caused by baryons. The large volume and high resolution of MTNG provide a unique opportunity to perform this comparison over a broad halo-mass range with a large number of well-resolved satellites, thereby ensuring accurate measurements of halo profiles.

Matching central objects between a dark-matter-only and a full hydrodynamic simulation (both run with the same initial conditions) is quite simple, since we expect the total mass and position of the haloes to be similar. Unfortunately, this is not the case for satellite galaxies. Once a central becomes a satellite, due to several processes, including baryonic effects, the pair of satellites will no longer share either a similar position or mass. To match satellite galaxies robustly, we track each individual object to a previous snapshot until it becomes a central and has a mass close to its peak halo mass during its evolutionary history. This process allows us to successfully match more than 97.5% of the satellite galaxies above the minimum stellar-mass limit we set for this study. We find that satellite galaxies in the hydrodynamic simulation are ∼3% closer to the centre of their halo than subhaloes in the dark-matter-only simulation, and also move up to ∼3% slower (with a small dependence on their halo mass). We show that this translates into a difference of up to 10% in galaxy clustering at scales of ∼0.1 h⁻¹ Mpc. While definitely significant, these differences are much lower than those we would find when selecting galaxies using a threshold based on different galaxy properties, such as stellar mass, luminosity in a given band, or star formation rate. We therefore conclude that, while baryonic effects significantly affect galaxy clustering and should be taken into account in precision cosmology, their impact is superseded by selection effects, which should be the priority when improving the modelling of galaxies in dark-matter-only simulations.

We use a suite of zoom-in hydrodynamic simulations, each implementing modified baryonic physics, to assess how such changes affect satellite properties. Although the limited number of massive haloes prevents robust statistical constraints, our results indicate that, even under extreme models, the main trends reported in this paper are unlikely to vary by more than a factor of two.

The results of this work quantify the impact of baryons on the satellite distribution not only for the SHAM approach, but more broadly for any galaxy population or galaxy formation model based on subhalo occupation. This includes, for example, semi-analytical models (Kauffmann et al. 1993; Cole et al. 1994; Somerville & Primack 1999; Bower et al. 2006; Lagos et al. 2008, 2018; Benson 2010, 2012; Jiang et al. 2014; Croton et al. 2016; Henriques et al. 2020), conditional luminosity functions (e.g. Yang et al. 2003; Cacciato et al. 2013), subhalo conditional abundance matching (SCAM; Zheng et al. 2009; Guo et al. 2016), and semi-empirical models (Moster et al. 2020), among others. Methods based on the halo occupation number, which typically populate galaxies using dark-matter particles within haloes rather than resolved subhaloes, will also be affected by baryonic effects, but in a manner different from that quantified in this work.

The outline of this paper is as follows. Section 2 presents the MTNG simulation, as well as its dark-matter-only counterpart, and describes the galaxy and subhalo samples used, along with the galaxy-matching algorithm. The impact of the baryonic effects on the positions and velocities of subhaloes is shown in Sections 3 and 4, respectively. Section 5 shows how these differences in positions and velocities propagate into the galaxy correlation function. We finalise by presenting our conclusions in Section 6. Appendix A describes the algorithm used to modify the positions and velocities of subhaloes in the dark-matter-only simulation so as to emulate the impact of baryons, while Appendix B examines how our main results change under alternative implementations of baryonic physics.

Unless otherwise stated, the standard units in this paper are h⁻¹ M_⊙ for masses and h⁻¹ Mpc for distances. All logarithm values are in base 10.

2. Numerical simulation and galaxy samples

In this section, we present the cosmological hydrodynamic simulation used in this work, the MillenniumTNG (MTNG), as well as its dark-matter-only counterpart (Section 2.1). Section 2.2 describes the algorithm employed to match objects between the two simulations. Finally, in Section 2.3, we outline the criteria used to select galaxies in each sample for comparison.

2.1. The MTNG simulations

To quantify the impact of baryons on satellite galaxies, we require a hydrodynamic simulation with a large volume to ensure accurate statistics, as well as high resolution to resolve several satellites per halo. One of the few simulations available today with these characteristics is the largest hydrodynamic run of the MillenniumTNG suite of simulations (Pakmor et al. 2023; Hernández-Aguayo et al. 2023; Barrera et al. 2023; Hadzhiyska et al. 2023a,b; Delgado et al. 2023; Kannan et al. 2023; Bose et al. 2023; Contreras et al. 2023b). This suite comprises tens of dark-matter-only and hydrodynamic simulations, along with light-cone galaxy catalogues generated using semi-analytical models. The project extends the two well-known Millennium (Springel et al. 2005) and IllustrisTNG (Nelson et al. 2018; Springel et al. 2018; Marinacci et al. 2018; Pillepich et al. 2018; Naiman et al. 2018; Nelson et al. 2019) projects in a direction that enables accurate studies of the galaxy–halo connection and of the impact of baryonic physics on clustering, particularly at much larger cosmological volumes than previously possible.

The largest hydrodynamic simulation of this suite (which we refer to simply as the MTNG simulation) has a periodic volume of (500 h⁻¹ Mpc)³. It contains 4320³ dark matter particles and an equal number of gas cells, implying an average gas-cell mass of 2.00 × 10^7 h⁻¹ M_⊙ and a mass resolution of 1.12 × 10^8 h⁻¹ M_⊙ for the dark matter particles. The simulation adopts a Planck Collaboration XIII (2016) cosmology¹, identical to that used in the IllustrisTNG project.

The initial conditions for MTNG were generated by fixing the amplitudes of the initial power modes to their expected rms values (Angulo & Pontzen 2016), which substantially reduces the effect of cosmic variance, at least for second-order statistics. The simulation was run using the moving-mesh code AREPO (Springel 2010), which includes radiative cooling, star formation in the gas, the growth of supermassive black holes, and the associated feedback from both supernovae and active galactic nuclei. Galaxies form as agglomerations of star particles whose properties can be directly measured from the simulation. The run comprises 265 snapshots, ranging from z = 63 (when the initial conditions were generated) to z = 0. Due to space limitations, only the group catalogues were stored, except for a limited number of key snapshots that also included particle data. This number of group and subhalo outputs is several times greater than in most previous large-scale simulations, providing a unique opportunity to follow the detailed evolution of satellite galaxies within haloes.

To study the impact of baryons on satellite subhaloes, we need to match the galaxies from MTNG to a simulation with similar characteristics but composed only of dark matter. For this purpose, we use another simulation from the MillenniumTNG suite with the same number of dark matter particles, the same cosmology, and the same initial conditions and volume. Both simulations were run with the moving-mesh code AREPO and have identical snapshot outputs at the same redshifts. We refer to this simulation as MTNG-DMO. Although we use a simulation as similar as possible to the original MTNG run to minimise noise in our comparison, we also tested matching our galaxy sample to a lower-resolution simulation (with 2³ fewer particles) that shares the same cosmology, initial conditions, and number of snapshots as MTNG, but was run with the Gadget-4 code (Springel et al. 2021). We found only negligible differences that do not affect our main results.

2.2. Matched catalogues

Because MTNG and MTNG-DMO share identical initial conditions (ICs), most haloes in one simulation can identify a “partner” object in the other with similar mass, position, and formation history, which originated from the same overdensity in the ICs. To match the host haloes (not subhaloes or galaxies) between the two simulations, we employ a kd-tree algorithm from SciPy (Virtanen et al. 2020) to search for neighbouring objects in a four-dimensional space defined by position and logarithmic halo mass. We set the weight of being off by one dex in halo mass to be equal to being off by 1 h⁻¹ Mpc. We test other configurations, only finding minor differences and a noticeable change in our results. To ensure robust matched catalogues, the procedure is performed bi-directionally: we first search for haloes in MTNG-DMO that best match those in MTNG, and then repeat the search in the opposite direction.

It is not possible, however, to apply this same matching algorithm directly to satellite galaxies or subhaloes. Due to both numerical and physical effects (such as baryonic processes), satellites gradually deviate from their counterparts in the other simulation, and dynamical friction cause their masses to decrease at different rates. To overcome this, we use the merger trees of both simulations to trace the main progenitor branch of all satellite galaxies in MTNG-DMO back to the snapshot where they were still central objects, with masses of approximately 80% of their peak subhalo mass (i.e. the maximum mass reached during their evolution). We then repeat the matching procedure described above for the host haloes at this specific snapshot. This process is applied to all satellite galaxies residing in haloes more massive than 10¹² h⁻¹ M_⊙ and with peak maximum circular velocities above 100 km s⁻¹ (corresponding to the lowest cut in our selection samples; see Section 2.3 for details). We adopt 80% of the peak subhalo mass because, at this stage, the objects are still massive enough for reliable identification by our algorithm, yet not so close to the haloes into which they will eventually merge that the matching becomes ambiguous. We tested alternative thresholds and found only minor differences in our results. This technique successfully matches more than 97.5% of all objects in our sample. The remaining ∼2.5% is most likely due to rare failures of the merger-tree algorithms, which can occasionally identify or follow objects differently in the two independent simulations.

2.3. Target galaxy samples

To study the effects of baryons on satellite subhaloes, we selected subhaloes directly from the dark-matter-only simulation and matched them to their counterparts in the hydrodynamic simulation. We chose to select subhaloes from the dark-matter-only run for simplicity, as it contains fewer objects (probably caused by the lower mass of the stellar mass particle, and that they are usually found in the inner cores of the subhaloes), and because this facilitates the characterisation of baryonic effects in the type of simulation most commonly used to build mock galaxy catalogues. We also performed the analysis selecting objects first from the hydrodynamical simulation and found results consistent with those presented here. The subhaloes were selected according to the peak value of their maximum circular velocity, $V_{\max }\equiv \sqrt{GM(<r)/r}$, attained during their evolution (V_peak). As shown by Chaves-Montero et al. (2016), subhaloes selected by V_peak exhibit clustering properties similar to those of galaxies selected by stellar mass, making it a standard choice for creating mocks in SHAM-like models. As we will see in the following sections, baryonic effects are stronger in the inner regions of haloes. Selecting galaxies by stellar mass or V_peak also preferentially targets objects located in these inner regions (in contrast to star-forming samples, which are mainly found in the outer regions of haloes; Orsi & Angulo 2018). Consequently, this selection is expected to be sensitive to baryonic effects, making V_peak an appropriate variable for defining our samples.

The subhalo samples were constructed at z = 0 for three different number densities: 0.0316, 0.01, and 0.00316 h³ Mpc⁻³. We adopt z = 0 as our reference redshift because baryonic effects on the matter distribution are expected to be strongest at late times (Ferlito et al. 2023); Schaller:2025. In addition, the satellite fraction is generally higher at low redshift (e.g. Contreras et al. 2013), which reduces stochastic noise in our results. A graphical representation of these samples is shown in Fig. 1, where each sample corresponds to the subhaloes located to the right of the intersection between the horizontal lines and the cumulative V_peak function. To facilitate comparison with a galaxy sample selected by stellar mass, the upper panel of the figure indicates the equivalent stellar mass for subhaloes with log(V_peak) = 2.00, 2.25, 2.50, and 2.75, assuming the stellar mass function of MTNG and a one-to-one mapping between stellar mass and V_peak in the dark-matter-only simulation (i.e. the most massive galaxy in the hydrodynamic simulation is assigned the largest V_peak value in the dark-matter-only simulation, the second most massive galaxy is assigned the second largest value, and so on, as in a basic SHAM approach).

Fig. 1. Cumulative distribution of Vpeak for subhaloes in the MTNG-DMO simulation. We define Vpeak as the maximum circular velocity (Vmax) reached by a subhalo during its lifetime. The upper axis shows the corresponding stellar mass of these haloes after matching them to their counterparts in the hydrodynamic MTNG simulation.

For the final part of this paper, we also select galaxies directly from the hydrodynamic MTNG simulation based on their stellar mass (defined as the total mass of stellar particles within twice the stellar half-mass radius), luminosity in the r and U bands (computed as the summed luminosities of all stellar particles in the subhalo), and star formation rate (SFR). These samples are defined to have the same number densities as the subhalo samples selected by V_peak.

3. Impact of baryons on the position of galaxies

We now examine the positional differences between the V_peak-selected satellite subhaloes of MTNG-DMO and their matched galaxies in the hydrodynamic MTNG simulation. The upper panel of Fig. 2 shows the mean distance of satellites from their halo centres as a function of halo mass for the n_den = 0.0316 h³ Mpc⁻³ sample. The solid line represents the median of the distribution, and the shaded region indicates the 16–84% range. Satellite galaxies in the hydrodynamic simulation are found to be 3–4% closer to their halo centres than their counterparts in the dark-matter-only simulation. A reduction in satellite distances is consistent with a scenario in which satellites lose kinetic energy and reduce their orbital angular momentum as they fall inward due to dynamical friction (Ogiya & Burkert 2016). Ram-pressure forces also contribute for galaxies that retain some of their gas, since gas experiences a stronger drag than dark matter. Through gravitational coupling within the satellite, this force is effectively transmitted to the system as a whole. In addition, satellites in the hydrodynamic simulation experience gaseous drag (i.e. the general retarding force on an object moving through gas), which differs from ordinary dynamical friction against the dark-matter background. For slightly supersonic motion – a regime commonly encountered by satellites – the gas drag can be even stronger (Ostriker 1999). Two other effects that can cause differences between a dark-matter-only run and one including baryons are tidal shocking – in particular, disk shocking, which can lead to satellite depletion in the inner regions of halos (D’Onghia et al. 2010) – and modifications to the host’s density profile and the associated strengthening of tidal forces due to the presence of a central galaxy.The combination of these mechanisms leads to a progressive loss of orbital angular momentum, which becomes most noticeable after each pericentric passage (Miller et al. 2020). Consistently, when examining individual orbits, we find that the discrepancy between the simulations emerges after the first pericentric passage. An example of three such satellite orbits is shown in Fig. 3.

Fig. 2. Top: median (solid lines) and 16th–84th percentile region (shaded area) of the mean radial distance between satellite galaxies and the halo centre as a function of halo mass. Results from the MTNG-DMO simulation are shown in black, and those from the hydrodynamic MTNG run in orange. Subhaloes were selected according to their Vpeak for a number density of nden = 0.0316 h3 Mpc−3 in MTNG-DMO, and the corresponding matched galaxies in MTNG were identified following the procedure described in Section 2.2. Distances are shown as a function of the MTNG-DMO halo masses. The black dashed line represents the value of one virial radius r200c. Bottom: Median (solid lines) and 16th–84th percentile region (shaded area) of the ratio of mean distances between matched haloes in the MTNG and MTNG-DMO simulations. The orange solid line and shaded area correspond to the galaxy sample with number density 0.0316 h3 Mpc−3 (same as in the top panel), while the dashed and dotted lines represent galaxy samples with number densities of 0.01 and 0.00316 h3 Mpc−3. Only the subhaloes successfully matched between the two simulations are shown here.

Fig. 3. Illustration of three satellite galaxies infalling into a halo of ∼ 1013.5 h−1 M⊙. Black circles show the satellite positions in the MTNG-DMO simulation, while orange circles indicate their matched counterparts in the hydrodynamic MTNG run. The symbol sizes scale with the mass of the satellites. The dashed grey line represents one virial radius.

Figure 2 also shows no significant dependence of this offset on halo mass. The bottom panel of the figure, which presents the ratio between the mean satellite distances in the MTNG and MTNG-DMO simulations, includes the lower-density samples and reveals no notable change in trend across different number densities. To minimise noise, we restrict the analysis in this panel to haloes that host, on average, at least five satellites per halo (i.e. ⟨N_sat(M_h, min)⟩ > 5). Lower-mass haloes display similar behaviour (not shown here).

We next explore how the baryonic effects depend on distance from the halo centre. Figure 4 shows the mean satellite density profiles for haloes with log(M_h/h⁻¹ M_⊙) ≃ 13.5, 14, and 14.5 in the MTNG and MTNG-DMO simulations. The halo mass bins have a width of 0.1 dex. As in the previous figure, the abundance of galaxies in the hydrodynamic simulation is higher at smaller radii than in the dark-matter-only run. The difference between simulations becomes most apparent on small scales, where the abundance of satellites in MTNG exceeds that in MTNG-DMO by up to ∼10% at radii of ∼30 h⁻¹ kpc, largely independent of halo mass. Only the highest-density sample is shown here, as the other samples exhibit similar trends. These results are consistent with those found by Vogelsberger et al. (2014b) and Vogelsberger et al. (2014a) in the original Illustris hydrodynamic simulation.

Fig. 4. Density of satellite galaxies as a function of distance from the halo centre for the nden = 0.0316 h3 Mpc−3 sample. The left, middle, and right panels show density profiles for halo masses of 1013.5, 1014, and 1014.5 h−1 M⊙, respectively. Black lines represent results from the MTNG-DMO simulation, while orange lines correspond to the hydrodynamic MTNG simulation. The vertical dashed lines represent the value of one virial radius r200c.

The enhanced abundance of satellites in the inner regions of host haloes may result from the increased central density – particularly of gas – on these scales. As satellites pass through pericentre, they lose orbital energy and slow down, bringing them closer to the halo centre and increasing the time they spend in this region, thereby affecting their radial distribution. We will examine this interpretation further in the next section.

4. Impact of baryons on the velocities of galaxies

We now focus on the velocities of satellite galaxies within haloes. To quantify the differences between simulations, we compute the relative velocity of each satellite with respect to its host halo, v′_sat = v_sat − v_h. We then define the satellite velocity dispersion as the three-dimensional velocity dispersion of the satellites within a halo.

Figure 5 shows the median and 16–84% range of the mean satellite velocity dispersion in MTNG and MTNG-DMO as a function of halo mass. In contrast to the behaviour observed for distances to the halo centre, the mean velocity dispersion is nearly identical for the most massive haloes. For less massive haloes, we find small differences of a few per cent, with satellites in the dark-matter-only simulation moving slightly faster than in the hydrodynamic run. The agreement across the three number densities (bottom panel) is even better than for positions, showing no systematic trend across samples.

Fig. 5. Top: median (solid lines) and 16th–84th percentile region (shaded area) of the satellite velocity dispersion as a function of halo mass for a number density of nden = 0.0316 h3 Mpc−3. Results from the MTNG-DMO and MTNG simulations are shown in black and orange, respectively. Bottom: median (solid lines) and 16th–84th percentile region (shaded area) of the ratio of velocity dispersions between matched haloes in the MTNG and MTNG-DMO simulations, for the three number densities used in this work (as labelled).

Although the overall velocity dispersion is very similar between the two simulations, a clear difference appears when considering the dependence on distance from the halo centre, as shown in Fig. 6. In the inner regions of haloes, satellite velocities are up to ∼10% lower in the hydrodynamic simulation. This difference is primarily caused by dynamical friction that slow down satellites as they pass through pericentre.

Fig. 6. Mean velocity dispersion of satellite galaxies as a function of distance from the halo centre for the nden = 0.0316 h3 Mpc−3 sample. The left, middle, and right panels show velocity profiles for halo masses of 1013.5, 1014, and 1014.5 h−1 M⊙, respectively. Black lines represent results from the MTNG-DMO simulation, while orange lines correspond to the hydrodynamic MTNG run. The vertical dashed lines represent the value of one virial radius r200c.

The decrease in velocity within the central regions is present at all halo masses, which may seem at first inconsistent with Fig. 5, where we find comparable mean velocity dispersions for the most massive haloes. This apparent discrepancy arises because satellites in the hydrodynamic simulation move more slowly in the inner regions but remain there for longer periods, as seen in Fig. 4. Although their instantaneous velocities are lower, these satellites still move faster than those in the outer halo regions. The two effects compensate, keeping the mean satellite velocity nearly constant for massive haloes and producing a small increase in velocity in the dark-matter-only simulation towards lower masses. This interpretation is consistent with feedback processes being more effective at lower halo masses – particularly supernova feedback, while AGN feedback remains active (Akino et al. 2022) – and thus with stronger baryonic effects in the inner parts of the mass density profile. After the first pericentric passage, satellites retain slightly lower velocities, as shown in the 0.1–1.0 h⁻¹ Mpc range of Fig. 6. In this regime, the population includes recently infalling haloes, so the velocity suppression is less pronounced. Beyond ∼1 h⁻¹ Mpc, the velocity distributions of the two simulations converge, as these objects are about to begin infall or are nearby systems identified as satellites by the halo finder, which we do not consider further in this analysis.

We conclude from these results that baryons directly affect satellite velocities, with the strongest influence occurring near pericentric passage. This effect also modifies the spatial distribution of satellites within haloes, drawing them closer to the centre and increasing the time they spend in the inner regions. Because satellites move faster in the inner parts of haloes, the median (and mean) velocity dispersion of the system does not decrease substantially, even though individual satellites experience lower orbital speeds along their trajectories. Therefore, this statistic should be interpreted with care when assessing the impact of baryons on satellite dynamics.

5. Baryonic and selection effects on galaxy clustering

In this section, we examine the impact of baryonic effects on several two-point statistics. Section 5.1 describes how we modified the dark-matter-only simulation to reproduce the baryonic effects observed in the hydrodynamic run and presents their impact on subhalo clustering. A more technical description of the baryon-mimicking procedure is provided in Appendix A. In Section 5.2, we compare the magnitude of baryonic effects with those introduced by galaxy-selection criteria.

5.1. The impact of baryons on galaxy clustering

In the previous sections, we explored the direct impact of baryons on satellite positions and velocities, finding variations of up to 10% on small scales. We now focus on how these differences affect galaxy clustering in real and redshift space. So far, our analysis has been based on subhaloes in MTNG-DMO above 10¹² h⁻¹ M_⊙ that could be matched to the hydrodynamic MTNG. Although the completeness of this matching exceeds 97.5%, computing clustering from an incomplete sample can introduce systematic biases. To account for this, we parameterise the positional and velocity differences found in Sections 3 and 4, and apply corresponding corrections to ensure that all subhaloes in the dark-matter-only simulation follow the same trends as in the hydrodynamic run. Below, we briefly describe these corrections; further details are given in Appendix A.

5.1.1. Change in the positions of satellites

We shift the positions of satellite subhaloes in MTNG-DMO to reproduce the radial distribution of galaxies in MTNG. The cumulative satellite density profile is computed for both simulations as a function of distance from the halo centre, normalised by the virial radius. This statistic is found to be nearly independent of halo mass. We then derive the mean radial displacement as a function of cumulative satellite abundance and apply this relation to all satellite subhaloes in MTNG-DMO, effectively reducing their distances from the halo centres.

5.1.2. Change in the velocities of satellites

To correct for velocity differences, we measure the ratio of satellite velocity dispersions between the simulations as a function of halo mass (analogous to Fig. 6). This ratio is well approximated by an error-function-like expression, equal to 1 on large scales and approaching 0.9 at small scales, with a smooth transition between the two. The small-scale value shows a weak dependence on halo mass, which we include in the model. We then apply this function to all subhaloes in the dark-matter-only simulation, thereby reducing their velocities to match the distribution in the hydrodynamic run. It is essential to apply this velocity correction after the positional shift; otherwise, the median velocity of the sample cannot be correctly recovered, as discussed in Section 4.

5.1.3. Change in the velocities of centrals

While the positions of central subhaloes and central galaxies coincide with the halo centre (i.e. the potential minimum), their velocities differ from that of the halo itself – typically defined as the mass-weighted mean velocity of all particles or cells in the group. This offset, known as the velocity bias, correlates with halo velocity dispersion and can significantly influence galaxy clustering (Guo et al. 2015). By construction, our use of central subhaloes from the dark-matter-only simulation implicitly includes velocity bias; however, its amplitude may differ from that in the hydrodynamic simulation. To create an equivalent galaxy sample and minimise systematic differences, we correct for this bias.

We first match all haloes above 10¹¹ h⁻¹ M_⊙ (the minimum halo mass in our densest sample) and measure the velocity difference as a function of halo mass. We define the velocity bias as the magnitude of the velocity difference between the host halo and its central subhalo, Δv_cen. We then compute the ratio of the mean Δv_cen between the two simulations and parameterise this ratio as a function of halo mass in MTNG-DMO, which allows us to correct all haloes accordingly.

The algorithm developed to modify the positions and velocities of central and satellite subhaloes can be applied to any dark-matter-only simulation to emulate the baryonic effects observed in MTNG. These modifications, however, do not constitute a universal prescription for baryonification; they merely reproduce the behaviour specific to this galaxy-formation model. Baryonic effects vary across simulations, and their amplitude in the real Universe remains uncertain (e.g. Aricò et al. 2020, 2021; see also Appendix B for a discussion of how different galaxy-formation models affect subhalo positions and velocities). Furthermore, although baryonification models typically depend on halo properties such as distance to the halo centre and halo mass (as we do here), secondary dependencies on halo concentration or large-scale environment may also play a role (Wang et al., in prep.). The purpose of this work is not to build a general baryonification model for subhaloes but rather to quantify, to first order, how baryons affect galaxy clustering. Nevertheless, the methodology presented here provides a useful framework for developing more universal subhalo baryonification schemes.

We refer to MTNG-DMO_mod as the MTNG-DMO simulation with subhalo positions and velocities modified according to the prescriptions above. In Fig. 7, we compare the clustering of subhaloes in MTNG-DMO and MTNG-DMO_mod. We compute the real-space two-point correlation function (ξ(r)), the projected correlation function in redshift space with an integration limit of 30 h⁻¹ Mpc (w_p(r)), the monopole of the redshift-space correlation function (ξ_ℓ = 0(r)), and the quadrupole (ξ_ℓ = 2(r)) for the three number-density samples used in this work. All clustering measurements were obtained with CORRFUNC (Sinha & Garrison 2019). To test the validity of this approach, we measure the clustering of the MTNG and the MTNG-DMO_mod using the matching catalogues, finding an overall good agreement in the clustering (see Appendix A for more details).

Fig. 7. Top: from left to right: the real-space two-point correlation function, the projected correlation function, the monopole, and the quadrupole of the correlation function for samples with number densities of 0.0316 (black lines), 0.01 (red lines), and 0.00316 h3 Mpc−3 (blue lines). Solid lines represent the standard MTNG-DMO samples, while dashed lines correspond to the clustering from the modified MTNG-DMO simulation, where satellite positions and velocities were adjusted to mimic those of the hydrodynamic MTNG run (MTNG-DMOmod). Bottom: ratio of the clustering between MTNG-DMOmod and MTNG-DMO for the two-point, projected, and monopole correlation functions. For the quadrupole (rightmost panel), we show the difference rather than the ratio to avoid discontinuities when the signal crosses zero.

The real- and redshift-space correlation functions, as well as the projected correlation function, show clustering differences of 10–20% on small scales, which gradually diminish and become negligible at scales of 1–3 h⁻¹ Mpc (bottom panels of Fig. 7). Chaves-Montero et al. (2016) also studied this effect for galaxies selected in bins of stellar mass and found compatible results. As expected from previous results, subhaloes in MTNG-DMO_mod are more strongly clustered than in MTNG-DMO. The impact of the velocity correction is modest for the monopole, producing differences similar to those seen in the real-space two-point function. Although clustering amplitudes vary among the number-density samples, the relative ratios – interpreted as the effect of baryons on galaxy clustering – show no significant dependence on number density.

The quadrupole of the correlation function, on the other hand, exhibits a clear dependence on number density: less dense samples (corresponding to more massive or luminous galaxies in an abundance-matching context) show larger amplitude differences than the densest samples. This behaviour, however, partly reflects the choice of presentation. Unlike the other statistics, where we plot ratios, here we show the difference in the signal multiplied by r², as the quadrupole crosses zero and the ratio would be ill-defined. When inspecting the ratio directly (not shown here), the number-density dependence weakens considerably. In any case, as discussed in the following section, the quadrupole amplitude is the least sensitive of the statistics considered here and does not warrant further interpretation at this stage. We therefore conclude that the effect of baryons on galaxy clustering is similar across all number densities examined in this work.

5.2. Baryonic effects and selection effects

Baryonic effects are not the only systematics that produce differences in the spatial distributions of satellite galaxies within haloes. Selection effects – i.e. the particular criteria used to define a galaxy sample in a survey – also modify the spatial and kinematic properties of satellite populations. For instance, the satellites in a star-forming sample tend to reside in the outer regions of haloes (Orsi & Angulo 2018), whereas galaxies selected by stellar mass are predominantly located in the inner regions. To account for both selection and baryonic effects, modern abundance-matching approaches (e.g. Contreras et al. 2021b; Favole et al. 2022; Rocher et al. 2023; Ortega-Martinez et al. 2024) and halo-occupation-based models (e.g. Carretero et al. 2015; Yuan et al. 2022; Gonzalez et al., in prep.) include parameters that explicitly control the spatial and/or velocity distribution of satellite galaxies. This can be implemented either by selecting dark-matter particles or subhaloes in specific regions of the halo (inner or outer) and populating them with galaxies, or by assigning galaxies directly through analytical prescriptions.

The dependence of galaxy clustering on stellar mass, different magnitude bands, and star-formation rate is well established in the literature. More massive or brighter galaxies, which typically populate more massive haloes, are more strongly clustered than less massive systems (e.g. Zehavi et al. 2011; Guo et al. 2013). The small-scale clustering of red galaxies is higher than that of blue galaxies, reflecting their different satellite fractions and radial distributions within haloes (e.g. Zehavi et al. 2005; Coil et al. 2006; Hartley et al. 2013). A similar behaviour is found when selecting galaxies by specific star-formation rate: systems with lower values of this quantity exhibit stronger small-scale clustering than highly star-forming galaxies (e.g. Li et al. 2006; Sobral et al. 2010). These empirical trends motivate our comparison of clustering statistics for samples selected in different ways, and provide observational context for the baryonic effects identified in cosmological simulations.

To place the impact of selection effects in context with that of baryonic effects, we compare the ratios and differences of several two-point statistics: the real-space and projected correlation functions, and the monopole and quadrupole of the redshift-space correlation function. These quantities are evaluated for galaxies selected by stellar mass (M_stell > 2.6 × 10⁸ h⁻¹ M_⊙), and for those selected by their r-band magnitude (M_r < −16.98), U-band magnitude (M_U < 16.75), and star formation rate (SFR > 0.0166 M_⊙/yr). We define the stellar mass as the total mass of all stellar particles within the stellar half-mass radius, the SFR as the sum of the individual star-formation rates of all gas cells associated with each subhalo, and the magnitudes from the summed luminosities of all stellar particles belonging to the group. All samples are drawn from the hydrodynamical MTNG simulation and correspond to galaxies with the highest stellar masses and SFRs, and the lowest magnitudes (i.e. highest luminosities), at a fixed number density of 0.0316 h³ Mpc⁻³ (Fig. 8), in direct analogy to the V_peak selection applied to the subhalo catalogue of the MTNG-DMO simulation.

Fig. 8. Ratio and difference in galaxy clustering between galaxies selected by Mr luminosity (dark red lines), Mu luminosity (dark blue lines), and star formation rate (SFR, magenta lines), and those selected by stellar mass in the MTNG simulation. Each galaxy sample has a number density of 0.0316 h3 Mpc−3. For comparison, dashed grey lines show the ratio and difference in clustering between Vpeak-selected subhaloes in the MTNG-DMOmod and MTNG-DMO simulations (ϕ/ϕmod, equivalent to the black solid line in the bottom panel of Fig. 7).

We find that the clustering differences among the selection-based samples are substantially larger than those induced by baryonic effects. On large scales, the projected and real-space correlation functions – both with and without redshift-space distortions – exhibit constant relative offsets. These reflect differences in the large-scale bias of the galaxy samples, which arise from variations in the satellite fraction and in the mass distribution of their host haloes (the dominant effect), as well as from different levels of galaxy assembly bias (Gao & White 2007; Croton et al. 2007; Zehavi et al. 2018; Artale et al. 2018; Contreras et al. 2019; García-Moreno & Chaves-Montero 2026).

For the real-space and projected correlation functions, we observe strong differences in the one-halo term, which vary in the opposite direction to the shifts produced by baryonic effects. As previously discussed, moving from redder samples (stellar-mass-selected galaxies) to bluer samples (M_r, M_U, and SFR) selects galaxies that increasingly occupy the outer regions of haloes, thereby reducing the small-scale clustering amplitude. This trend is not present for the monopole of the correlation function. We interpret this as evidence that selecting outer-halo galaxies also selects systems with different velocity profiles: satellites in the outskirts have lower orbital velocities and follow more coherent orbits around the halo centre. As a result, they appear dynamically closer to the centre than the faster-moving satellites in the inner regions. These velocity-profile differences become more evident in the quadrupole, which shows larger deviations as we move towards bluer samples. The competition between spatial and velocity effects reduces the overall difference in the monopole at the smallest scales.

The differences we find among the various selection functions are far greater than those caused by baryonic effects (indicated by the grey dashed lines in the figures). This is particularly evident for the quadrupole, where the selection-induced deviations exceed the baryonic ones by more than an order of magnitude. We also tested different number densities and enforced equal satellite fractions across samples (not shown), finding some variation in the detailed trends but reaching the same conclusion: selection effects have a significantly stronger impact on galaxy clustering than baryonic effects.

6. Summary

In this work, we have examined the impact of baryonic effects on the positions and velocities of satellite galaxies. We used galaxies from MTNG, one of the largest high-resolution cosmological hydrodynamic simulations available, and matched them to subhaloes from a dark-matter-only simulation, MTNG-DMO, which shares the same initial conditions, cosmology, and resolution as its hydrodynamic counterpart. Satellites in both simulations were matched by tracing their progenitors to the epoch when they were central objects and selecting those with similar positions and halo masses. This procedure ensures a reliable correspondence between simulations, even for satellites that fell into their host haloes several gigayears ago and whose present-day positions and masses differ substantially.

Here we summarise the main results of this work:

• Satellite galaxies in the hydrodynamic simulation are, on average, 3–4% closer to the halo centre than the subhaloes in the dark-matter-only simulation (Fig. 2). When analysing the radial profile of satellite distribution, most of the difference originates from objects located in the innermost halo regions (100–300 h⁻¹ kpc), where the hydrodynamic simulation exhibits up to a 10% higher abundance (Fig. 4). We find no significant dependence on halo mass.

• The median satellite velocity dispersion of haloes above 10¹⁴ h⁻¹ M_⊙ is similar in both simulations. For lower-mass haloes, satellites in the dark-matter-only simulation move up to 3% faster than in the hydrodynamic run (Fig. 5). When considering the dependence of velocity on distance from the halo centre, galaxies within 100–200 h⁻¹ kpc show the largest differences, with satellites in the dark-matter-only simulation moving up to 10% faster than those in the hydrodynamic simulation. These differences show little correlation with halo mass (Fig. 6). We do not observe a change in the median velocity dispersion for the most massive haloes in the hydrodynamic run, since although satellites move more slowly in the inner regions, their higher abundance there compensates for this effect. Because satellites move faster in the inner halo, the combination of slower velocities and higher density yields a similar overall velocity dispersion.

• We developed an analytic approach to modify the positions and velocities of subhaloes in a dark-matter-only simulation so that they mimic the baryonic effects observed in MTNG (Section 5 and Appendix A). The corrections are parameterised as functions of the satellite distance from the halo centre and host-halo mass.

• We compared the subhalo clustering in MTNG-DMO with that in MTNG-DMO_mod, a modified version of the dark-matter-only simulation in which subhalo positions and velocities were adjusted to reproduce the baryonic effects seen in the hydrodynamic MTNG run. We did not directly compare MTNG and MTNG-DMO because the matching between the two is not fully complete, and differences in selection functions could introduce systematic biases, but we test the validity of our approach by comparing the clustering of the matched catalogues and find good agreement with the MTNG. We find that MTNG-DMO_mod exhibits higher clustering amplitudes in the two-point, projected, and monopole correlation functions – by up to 10–20% at scales of ∼100 h⁻¹ kpc – with the differences vanishing at scales of ∼1 h⁻¹ Mpc. We find no significant dependence on number density for V_peak-selected subhalo samples. The quadrupole shows only minor differences between models, with Δr²ξ_ℓ = 2 ≈ 2 h⁻² Mpc² (Fig. 7).

• We also compared the clustering of galaxies in MTNG selected by different intrinsic properties, including stellar mass, M_r, M_U, and SFR (Fig. 8). The clustering differences among these selections are much larger than those caused by baryonic effects. Only the monopole shows differences of comparable magnitude at small scales, likely due to compensating effects between galaxy positions and velocities. More complex colour-based selection functions, such as those used in surveys like DESI or Euclid (e.g. Raichoor et al. 2023; Myers et al. 2023; Euclid Collaboration: Lesci et al. 2024), may produce even larger differences.

We emphasise that the baryonic effects quantified here correspond to the fiducial MTNG galaxy-formation model and may not precisely reflect those in the real Universe. In fact, the baryonic suppression of the mass distribution in MTNG is somewhat weaker than that predicted by other cosmological hydrodynamic simulations (e.g. Aricò et al. 2020, 2021). In Appendix B, we use a suite of 26 zoom-in hydrodynamic simulations with the same initial conditions as MTNG, each incorporating modified baryonic physics, to assess how such changes affect satellite properties. This set includes 452 re-simulated haloes (29 with M_h > 10¹³ h⁻¹ M_⊙). Different physical prescriptions yield different amplitudes for baryonic effects, with the most extreme models predicting an increase by up to a factor of two relative to the fiducial case. However, the limited number of massive haloes prevents robust statistical constraints; these results should thus be regarded as indicative, providing an upper bound on the expected amplitude of baryonic effects.

Even if baryonic effects are enhanced by a factor of two or three due to differences in galaxy-formation physics, they remain subdominant compared to selection effects. Nevertheless, they are not negligible, as they correspond to 1–3σ differences in clustering assuming DESI 1% errors (Yuan et al. 2024). These discrepancies will be larger for DR1 or other future releases of DESI. Therefore, baryonic effects should be explicitly accounted for in galaxy-modelling frameworks.

Overall, our results provide a detailed view of how baryons influence position and velocities of satellite galaxies and how these effects depend on halo mass and radial position. These findings will contribute to improving how baryonic corrections are implemented in empirical models, particularly in conjunction with selection effects. Although developing a universal baryonification model was not the goal of this paper, the method introduced here offers a promising basis for extending baryonic corrections to any dark-matter-only simulation. In future work, we will investigate in greater detail how these effects depend on galaxy-formation physics and test the performance of empirical models in reproducing them.

Acknowledgments

We thank Giovanni Arico and Jonas Chaves-Montero for some usfull suggestions and comments. SC acknowledges the support of the “Ramón y Cajal” fellowship (RYC2023-043783-I). REA received support from grant PID2024-161003NB-I00 funded by MICIU/AEI/10.13039/501100011033 and by ERDF/EU. SB is supported by the UKRI Future Leaders Fellowship [grant numbers MR/V023381/1 and UKRI2044]. FM acknowledges support by the Simons Collaboration on “Learning the Universe”. VS acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. Author Contributions Statement: The idea for this project was developed by SC, with essential contributions from REA. SC analysed the data, performed the calculations, and wrote the manuscript. REA provided critical comments throughout the development of the project. SB, BH, LH, RP, and VS developed the MTNG simulation. RP provided critical feedback on Appendix A, and VS offered technical assistance during the analysis of the simulations. REA, FM, and VS developed the zoom-in simulations used in Appendix B, and FM ran most of them. All authors contributed comments on the manuscript. From the third author onward, the author list is ordered alphabetically.

References

Appendix A: Mimicking baryonic effects in a dark-matter-only simulation

In this appendix, we detail the algorithm used to modify the positions and velocities of subhaloes in MTNG-DMO so as to reproduce the baryonic effects present in the hydrodynamic MTNG run. The procedure is divided into three steps. We first modify the positions of satellite subhaloes, then adjust their velocities, and finally correct the velocities of central galaxies.

A.1. Modifying the satellite positions

To modify satellite positions, we first compute the cumulative number of satellites as a function of their distance to the halo centre, normalised by the virial radius r_200c. We normalise this distribution by the total number of satellites per halo, which restricts the function to the interval [0, 1]:

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{c}}(d_{\mathrm{radial}}/r_{200c})=\frac{N_{\mathrm{sat}}(<d_{\mathrm{radial}}/r_{200c})}{N_{\mathrm{sat},\mathrm{tot}}},\end{array}$$(A.1)

where N_sat(< d_radial/r_200c) is the number of satellites within d_radial and N_sat, tot is the total number of satellites in the halo.

Normalising distances by r_200c and the counts by N_sat, tot yields a distribution with very weak halo-mass dependence. We then compute the difference between the inverse cumulative functions of the hydrodynamic and dark-matter-only simulations, i.e. the difference in the average satellite distance at fixed cumulative satellite fraction:

$$\begin{array}{c}\hfill \mathrm{\Delta }{r}_{\mathrm{sat}}(n_{\mathrm{sat},\mathrm{rad}})={\varphi }_{c,\mathrm{hydro}}^I(n_{\mathrm{sat},\mathrm{rad}})-{\varphi }_{c,\mathrm{dmo}}^I(n_{\mathrm{sat},\mathrm{rad}}),\end{array}$$(A.2)

where Δr_sat(n_sat, rad) is expressed in units of the virial radius and

$$\begin{array}{c}\hfill n_{\mathrm{sat},\mathrm{rad}}=\frac{N_{\mathrm{sat}}(<d_{\mathrm{radial}}/r_{200c})}{N_{\mathrm{sat},\mathrm{tot}}}\end{array}$$(A.3)

is the satellite rank within the halo (i.e. n_sat, rad = 1 for the outermost satellite and n_sat, rad = 1/N_sat, tot for the innermost satellite).

We find that Δr_sat(n_sat, rad) shows almost no dependence on host-halo mass. We therefore compute an average Δr_sat(n_sat, rad) combining all haloes with M_h > 10¹³ h⁻¹ M_⊙. This is done using only the matched subhaloes (about 97% of the total sample for our highest number-density case). Finally, we modify the satellite distances as

$$\begin{array}{c}\hfill d_{\mathrm{radial}}^{\prime }=d_{\mathrm{radial}}+\mathrm{\Delta }{r}_{\mathrm{sat}}(n_{\mathrm{sat},\mathrm{rad}}).\end{array}$$(A.4)

The correction is applied only to the modulus of the distance to the halo centre; we do not change the direction of the position vector, nor do we alter the spatial alignment of satellites.

We calibrated Δr_sat(n_sat, rad) using the sample with number density 0.0316  h³ Mpc⁻³, which shows the lowest noise. The resulting satellite distributions for MTNG-DMO_mod, MTNG, and MTNG-DMO are shown in the upper panel of Fig. A.1. This correction reproduces well the radial distributions over the full halo-mass range. We also checked the median satellite distance as a function of halo mass (analogous to Fig. 2) and found a similar level of accuracy (not shown).

Fig. A.1. Similar to Figs. 4 and 6, but including the predictions of the MTNG-DMOmod simulation (dashed lines), in which the positions and velocities of subhaloes were modified to mimic those of the hydrodynamic MTNG simulation.

A.2. Modifying the satellite velocities

We now turn to satellite velocities. We start by measuring, in both simulations, the mean modulus of the satellite velocity relative to the host halo, as a function of the satellite distance to the halo centre and in bins of halo mass. At this point, we use the satellite positions already corrected as described above. It is essential to correct positions first, because satellite velocity strongly correlates with radius; if the order is inverted, the median velocity of the final sample will not match that of the hydrodynamic run.

We denote the average satellite speed in a given simulation as

$$|v|(d_{\mathrm{radial}},M_{\mathrm{h}})\equiv ⟨|{\mathbf{v}}_{\mathrm{sat}}-{\mathbf{v}}_{\mathrm{h}}|⟩,$$

and consider the ratio between the hydrodynamic and dark-matter-only simulations. We find that, at small radii, this ratio is ∼0.9, while at large radii it approaches 1, with a smooth transition between both regimes. To model this behaviour and apply it to all subhaloes, we describe the ratio with an error-function-like parametrisation:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{|v_{\mathrm{hydro}}|}{|v_{\mathrm{dmo}}|}(d_{\mathrm{radial}})=& \frac{(f_{max}-f_{min})}{2}\mathrm{erf}\left(\frac{log(d_{\mathrm{radial}})-log(d_0)}{\sigma }\right)\hfill \\ \hfill & +\frac{(f_{max}+f_{min})}{2},\hfill \end{array}\end{array}$$(A.5)

where f_min and f_max are the minimum and maximum values of the ratio, d₀ is the characteristic radius of the transition, σ controls its smoothness, and erf is the error function,

$$\begin{array}{c}\hfill \mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x}{e}^{-t^2}\mathrm{d}t.\end{array}$$(A.6)

We find that f_min and f_max show a weak dependence on halo mass, which we model as

$$\begin{array}{c}\hfill f_{min}(M_{\mathrm{h}})=a_{f,min}(log{M}_{\mathrm{h}}-14)+f_{min,14},\end{array}$$(A.7)

and

$$\begin{array}{c}\hfill f_{max}(M_{\mathrm{h}})=a_{f,max}(log{M}_{\mathrm{h}}-14)+f_{max,14},\end{array}$$(A.8)

where f_min, 14 and f_max, 14 are the values for haloes of 10¹⁴ h⁻¹ M_⊙, and a_f, min and a_f, max encode the mild mass dependence.

For the MTNG hydrodynamic simulation we obtain the following best-fitting parameters:

$$(f_{min,14},f_{max,14},\sigma ,log{d}_0,a_{f,min},a_{f,max})=$$

$$(0.8733,1.004,0.285,-0.715,0.05,0.03).$$

The calibration was done using the highest number-density subhalo sample, which also agrees well with the lower-density samples.

The lower panel of Fig. A.1 shows the velocity-dispersion profiles for three halo-mass bins, comparing MTNG, MTNG-DMO, and MTNG-DMO_mod (grey dashed line). The corrected catalogue reproduces the hydrodynamic velocity profiles very well. We also tested a simplified correction that ignores the halo-mass dependence; this has negligible impact on both the velocity profiles and the clustering results of Fig. 7.

A.3. Change in the velocities of centrals

The subhalo finder SUBFIND assigns the same position to the central subhalo/galaxy and to the host halo (the potential minimum). Velocities, however, are usually defined as the mass-weighted mean of all particles or cells in the halo, so the velocity of the central object differs from that of the halo. This offset, known as velocity bias, scales with the halo velocity dispersion and can noticeably affect subhalo/galaxy clustering (Guo et al. 2015). Although a detailed study of baryonic effects on centrals is beyond the scope of this paper, we must correct for this velocity difference to isolate the impact of baryons on satellite clustering.

To do so, we extend our halo matching down to 10¹¹ h⁻¹ M_⊙, the minimum halo mass of our densest sample. Matching centrals between simulations with identical initial conditions is straightforward. We then compute, for each matched halo, the relative velocity of the central object with respect to the halo (central minus halo velocity) in both simulations and take the ratio. Figure A.2 shows this ratio as a function of halo mass. The largest differences occur around 10¹² h⁻¹ M_⊙ and 10^14.5 h⁻¹ M_⊙, while intermediate masses show smaller offsets.

Fig. A.2. Ratio of the relative velocities of central subhaloes with respect to their host haloes between the MTNG and MTNG-DMO simulations, as a function of host-halo mass. The grey dashed line shows a two-Gaussian parametrisation used to describe this relation and to correct dark-matter-only simulations.

We parametrise this behaviour as

$$\begin{array}{c}\hfill g(M_{\mathrm{h}})=1-G_1(log{M}_{\mathrm{h}})-G_2(log{M}_{\mathrm{h}}),\end{array}$$(A.9)

where G₁ and G₂ are Gaussian functions,

$$\begin{array}{c}\hfill G(x;a,b,c)=aexp[-\frac{{(x-b)}^2}{2c^2}].\end{array}$$(A.10)

The best-fitting parameters are

$$(a_1,b_1,c_1)=(0.139,11.93,0.376),$$

$$(a_2,b_2,c_2)=(0.166,14.48,0.732).$$

As shown in Fig. A.2, this parametrisation reproduces the velocity ratios over the full halo-mass range. We apply g(M_h) as a multiplicative correction to the central velocities of all haloes in the dark-matter-only simulation.

All the corrections above should be regarded as first-order. We did not explore possible dependencies on secondary halo properties (concentration, formation time, or environment), nor did we test whether the functional form persists for other galaxy-formation prescriptions. Nonetheless, the methodology can be used as a basis for a more accurate and universal implementation of baryonic effects in subhaloes.

As a proof-of-concept, we apply all corrections to the matched catalogue of the MTNG-DMO simulation and compare its clustering to that of the hydrodynamic MTNG run (Fig. A.3). We find very good agreement between the modified MTNG-DMO and the full MTNG simulation, with only minor differences at small scales in the multipoles of the correlation function. Although developing a full baryonification model for subhaloes is beyond the scope of this work, our results demonstrate that the relations presented here provide a solid foundation for constructing a more general and physically motivated subhalo baryonification framework.

Fig. A.3. Clustering measurements for the MTNG, MTNG-DMO, and the modified MTNG-DMO version (labelled MTNG-DMOmod). All results use the matched catalogues selected by a Vpeak threshold corresponding to a number density of 0.0316  h3 Mpc−3.

Appendix B: Dependence on the galaxy-formation physics

In this paper, we used the MTNG hydrodynamic simulation to quantify the effects of baryons on satellite positions and velocities and on galaxy clustering. While the trends we found may qualitatively hold for other galaxy-formation models, the amplitude of baryonic effects in the real Universe is not well constrained, and it is not obvious how baryonic effects on the matter distribution translate into baryonic effects on satellites. To explore the possible range of amplitudes, we analysed a suite of 26 zoom-in simulations that share the same initial conditions as MTNG but adopt variation of the galaxy-formation prescriptions. The initial conditions were produced using a multi-zoom technique described in Burger et al. (2025). Each of these simulations re-simulates 452 haloes at the same resolution as MTNG.

The zoom-in simulations were run with AREPO and adopt the same numerical setup, but differ in seven astrophysical parameters varied over broad ranges (Table B.1). The parameters were chosen using a Latin-hypercube sampling of the parameter space, such that all parameters were varied simultaneously in a way that maximizes coverage of the multidimensional hypervolume. Consequently, each of the 26 simulations explores a different combination of parameters. These parameters span feedback strengths, wind energetics, and star-formation efficiencies, and were chosen to cover a plausible range of galaxy-formation physics. This ensemble represents the first stage of a larger zoom-in programme (Maion et al., in prep.) and can be used to estimate an upper limit on baryonic effects in satellites.

We match haloes and subhaloes in each zoom-in run to those in the corresponding dark-matter-only zoom-in simulation following the same procedure as in Section 2.2. Each simulation contains 29 haloes above 10¹³ h⁻¹ M_⊙, which limits the statistical power of the analysis. To increase the signal, we consider all subhaloes with V_peak ≳ 76 km s⁻¹, corresponding to a number density of 0.1 h³ Mpc⁻³ in the full MTNG-DMO box.

Figure B.1 shows the ratio of the median satellite distance (top) and the ratio of the satellite velocity dispersion (bottom) between the zoom-in hydrodynamic runs and their dark-matter-only counterparts, analogous to Figs. 2 and 5. The zoom-in run with the fiducial MTNG physics is shown in black. We do not find strong trends with most galaxy-formation parameters, although there are weak correlations with parameters linked to galactic winds and to the star-formation timescale. In Fig. B.1, lighter (darker) orange lines indicate runs with lower (higher) values of WindEnergyIn1e51erg. There is a mild tendency for lower wind-energy models to place satellites closer to halo centres and to reduce their velocities consistent with expectations from enhanced orbit decay with increased baryonic density in the CGM, but the scatter is too large to draw firm conclusions. We remind the reader that these simulations vary all the astrophysical parameters simultaneously, which severely limits the ability to identify trends with a single parameter. Nevertheless, this study is useful to guide the design of more focused follow-up simulations that isolate the parameters showing the strongest impact on satellite properties.

The small number of massive haloes also limits our ability to constrain the maximum possible impact of baryons on satellites. We do see that baryonic effects can be up to a factor of ∼2 larger than in the fiducial model. Although part of this may be due to noise, this factor provides a reasonable first-order upper bound on how much baryons could shift satellite positions or velocities in other galaxy-formation scenarios. Because of the limited sample, we did not attempt to measure more demanding statistics. In future work, we will explore these dependencies using larger suites of hydrodynamic simulations.

Fig. B.1. Top: ratio of the median satellite distance to the halo centre as a function of halo mass between a suite of zoom-in hydrodynamic simulations and their dark-matter-only counterparts. The dark-matter-only haloes are selected with a Vpeak cut equivalent to a number density of nden = 0.1  h3 Mpc−3 in the full MTNG-DMO box. The zoom-in with the same physical model as MTNG is shown as a black dashed line. Orange lines show runs with varied galaxy-formation parameters; lighter (darker) colours correspond to lower (higher) values of WindEnergyIn1e51erg, the parameter that correlates best with the amplitude of the effect. Bottom: same as the top panel, but for the satellite velocity dispersion.

All Tables

All Figures

	Fig. 1. Cumulative distribution of Vpeak for subhaloes in the MTNG-DMO simulation. We define Vpeak as the maximum circular velocity (Vmax) reached by a subhalo during its lifetime. The upper axis shows the corresponding stellar mass of these haloes after matching them to their counterparts in the hydrodynamic MTNG simulation.
In the text

Fig. 2. Top: median (solid lines) and 16th–84th percentile region (shaded area) of the mean radial distance between satellite galaxies and the halo centre as a function of halo mass. Results from the MTNG-DMO simulation are shown in black, and those from the hydrodynamic MTNG run in orange. Subhaloes were selected according to their Vpeak for a number density of nden = 0.0316 h3 Mpc−3 in MTNG-DMO, and the corresponding matched galaxies in MTNG were identified following the procedure described in Section 2.2. Distances are shown as a function of the MTNG-DMO halo masses. The black dashed line represents the value of one virial radius r200c. Bottom: Median (solid lines) and 16th–84th percentile region (shaded area) of the ratio of mean distances between matched haloes in the MTNG and MTNG-DMO simulations. The orange solid line and shaded area correspond to the galaxy sample with number density 0.0316 h3 Mpc−3 (same as in the top panel), while the dashed and dotted lines represent galaxy samples with number densities of 0.01 and 0.00316 h3 Mpc−3. Only the subhaloes successfully matched between the two simulations are shown here.

In the text

	Fig. 3. Illustration of three satellite galaxies infalling into a halo of ∼ 1013.5 h−1 M⊙. Black circles show the satellite positions in the MTNG-DMO simulation, while orange circles indicate their matched counterparts in the hydrodynamic MTNG run. The symbol sizes scale with the mass of the satellites. The dashed grey line represents one virial radius.
In the text

Fig. 4. Density of satellite galaxies as a function of distance from the halo centre for the nden = 0.0316 h3 Mpc−3 sample. The left, middle, and right panels show density profiles for halo masses of 1013.5, 1014, and 1014.5 h−1 M⊙, respectively. Black lines represent results from the MTNG-DMO simulation, while orange lines correspond to the hydrodynamic MTNG simulation. The vertical dashed lines represent the value of one virial radius r200c.

In the text

Fig. 5. Top: median (solid lines) and 16th–84th percentile region (shaded area) of the satellite velocity dispersion as a function of halo mass for a number density of nden = 0.0316 h3 Mpc−3. Results from the MTNG-DMO and MTNG simulations are shown in black and orange, respectively. Bottom: median (solid lines) and 16th–84th percentile region (shaded area) of the ratio of velocity dispersions between matched haloes in the MTNG and MTNG-DMO simulations, for the three number densities used in this work (as labelled).

In the text

Fig. 6. Mean velocity dispersion of satellite galaxies as a function of distance from the halo centre for the nden = 0.0316 h3 Mpc−3 sample. The left, middle, and right panels show velocity profiles for halo masses of 1013.5, 1014, and 1014.5 h−1 M⊙, respectively. Black lines represent results from the MTNG-DMO simulation, while orange lines correspond to the hydrodynamic MTNG run. The vertical dashed lines represent the value of one virial radius r200c.

In the text

Fig. 7. Top: from left to right: the real-space two-point correlation function, the projected correlation function, the monopole, and the quadrupole of the correlation function for samples with number densities of 0.0316 (black lines), 0.01 (red lines), and 0.00316 h3 Mpc−3 (blue lines). Solid lines represent the standard MTNG-DMO samples, while dashed lines correspond to the clustering from the modified MTNG-DMO simulation, where satellite positions and velocities were adjusted to mimic those of the hydrodynamic MTNG run (MTNG-DMOmod). Bottom: ratio of the clustering between MTNG-DMOmod and MTNG-DMO for the two-point, projected, and monopole correlation functions. For the quadrupole (rightmost panel), we show the difference rather than the ratio to avoid discontinuities when the signal crosses zero.

In the text

Fig. 8. Ratio and difference in galaxy clustering between galaxies selected by Mr luminosity (dark red lines), Mu luminosity (dark blue lines), and star formation rate (SFR, magenta lines), and those selected by stellar mass in the MTNG simulation. Each galaxy sample has a number density of 0.0316 h3 Mpc−3. For comparison, dashed grey lines show the ratio and difference in clustering between Vpeak-selected subhaloes in the MTNG-DMOmod and MTNG-DMO simulations (ϕ/ϕmod, equivalent to the black solid line in the bottom panel of Fig. 7).

In the text

	Fig. A.1. Similar to Figs. 4 and 6, but including the predictions of the MTNG-DMOmod simulation (dashed lines), in which the positions and velocities of subhaloes were modified to mimic those of the hydrodynamic MTNG simulation.
In the text

	Fig. A.2. Ratio of the relative velocities of central subhaloes with respect to their host haloes between the MTNG and MTNG-DMO simulations, as a function of host-halo mass. The grey dashed line shows a two-Gaussian parametrisation used to describe this relation and to correct dark-matter-only simulations.
In the text

	Fig. A.3. Clustering measurements for the MTNG, MTNG-DMO, and the modified MTNG-DMO version (labelled MTNG-DMOmod). All results use the matched catalogues selected by a Vpeak threshold corresponding to a number density of 0.0316  h3 Mpc−3.
In the text

Fig. B.1. Top: ratio of the median satellite distance to the halo centre as a function of halo mass between a suite of zoom-in hydrodynamic simulations and their dark-matter-only counterparts. The dark-matter-only haloes are selected with a Vpeak cut equivalent to a number density of nden = 0.1  h3 Mpc−3 in the full MTNG-DMO box. The zoom-in with the same physical model as MTNG is shown as a black dashed line. Orange lines show runs with varied galaxy-formation parameters; lighter (darker) colours correspond to lower (higher) values of WindEnergyIn1e51erg, the parameter that correlates best with the amplitude of the effect. Bottom: same as the top panel, but for the satellite velocity dispersion.

In the text