© The Authors 2026

1. Introduction

The standard Λ cold dark matter (ΛCDM) cosmological model predicts that dark matter haloes grow hierarchically, with larger haloes assembling largely through mergers with smaller haloes. Following a halo merger, the galaxies originally associated with the smaller haloes eventually merge with the central galaxy of the larger halo (i.e. the galaxy at the centre of the halo, typically the most massive and brightest galaxy in the system) or be disrupted, although these events do not immediately occur. This scenario implies that massive haloes host multiple galaxies: one central galaxy and several satellite galaxies. During their infall into the host halo and over successive orbits, satellite galaxies are affected by several physical processes, including starburst episodes, tidal stripping, strangulation, and gas stripping, which significantly alter their properties.

The properties of satellite galaxies, as well as their abundance and distribution within their host haloes, are the result of the cosmological model and the physics of galaxy formation. Galaxies form and evolve within dark matter haloes (White & Rees 1978) until their host haloes merge with more massive systems, and the frequency of such halo mergers is set by the underlying cosmological model. Using N-body simulations, we have learned that the positions and velocities of satellite galaxies are well approximated by those of dark matter subhaloes in gravity-only simulations (i.e. dark matter-only simulations, whose evolution is determined solely by the cosmological model). This relation is one of the foundations of the galaxy–halo connection and has led to the development of the subhalo abundance matching (SHAM) family of models, which are widely used to construct mock galaxy catalogues (Vale & Ostriker 2006; Conroy et al. 2006; Chaves-Montero et al. 2016; Contreras et al. 2021b).

However, although the dynamics of galaxies are expected to be similar to those predicted by gravity-only simulations, they are not identical. In Contreras et al. (2026), we studied the impact of baryons on the positions and velocities of dark matter subhaloesusing the MillenniumTNG (MTNG) hydrodynamic simulation (Pakmor et al. 2023; Hernández-Aguayo et al. 2023), finding differences of a few per cent in the mean distance to the halo centre and in the velocity dispersion of satellite galaxies. Although this may appear to be a minor effect, Contreras et al. (2026) predict that it leads to 10–20% variations in the clustering amplitude, corresponding to 1–3σ for DESI-like errors on small scales (i.e. below 1 h⁻¹ Mpc). In addition, the detectability of satellite galaxies in observations depends on their baryonic content. As satellite galaxies lose their gas and undergo quenching, they can only be detected through observables linked to star formation if they were accreted recently, producing a significant detection bias towards the outskirts of their host dark matter haloes (Orsi & Angulo 2018).

A good understanding of how the properties of satellite galaxies evolve is crucial for improving empirical models based on the galaxy–halo connection, such as the halo occupation distribution (HOD; Jing et al. 1998; Benson et al. 2000; Peacock & Smith 2000; Berlind et al. 2003; Zheng et al. 2005, 2007; Contreras et al. 2013; Guo et al. 2015; Contreras et al. 2017; Contreras & Zehavi 2023) and SHAM. These models are widely used to characterise galaxy populations and to construct mock galaxy catalogues that play a central role in the analysis and interpretation of galaxy surveys (e.g. Rocher et al. 2023; Yuan et al. 2024; Ortega-Martinez et al. 2025), survey design and preparation (e.g. Smith et al. 2022), the construction of covariance matrices (e.g. Manera et al. 2013, 2015; Lin et al. 2020), and even the inference of cosmological parameters (e.g. Simha & Cole 2013; Contreras et al. 2023b; Mahony et al. 2026; Ortega-Martinez et al., in prep.). The way satellite galaxies are modelled in these mocks can significantly impact the resulting clustering predictions, particularly in the non-linear regime and, therefore, affect both cosmological and galaxy formation constraints derived from galaxy survey data.

In this work, we used the state-of-the-art hydrodynamic MTNG simulation to study the evolution of the baryonic content of satellite galaxies. In the first part of this paper, we focus on characterising the evolution of their stellar mass, gas mass, and r- and U-band magnitudes. At both z = 0 and z = 1, we describe the evolution of these satellite properties in terms of subhalo properties, including subhalo mass loss, pericentric distance (defined as the minimum distance between the satellite and the host-halo centre), time since infall (defined as the time elapsed since a galaxy was classified as a satellite), and the remaining maximum circular velocity fraction. Some of these subhalo properties are directly available in the gravity-only simulations commonly used to construct mock galaxy catalogues (e.g. Angulo et al. 2021), while others can be derived from their merger trees. We parametrised several of these relations to provide empirical models with a more accurate description of satellite galaxies. Some of these properties, particularly subhalo mass loss, have previously been linked to the evolution of satellite galaxies (e.g. Peñarrubia et al. 2008; Smith et al. 2013, 2016; He et al. 2026) and have also been used in empirical models of satellite evolution (e.g. Moster et al. 2020, who used subhalo mass loss, and Behroozi et al. 2019, who used changes in v_max). Nonetheless, a systematic study that explores the evolution of several galaxy properties as a function of multiple subhalo properties, while also assessing the performance and limitations of these approaches, remains absent from the literature.

In the second part of this paper, we investigate the evolution of the mass profiles of satellite galaxies during their infall into the host halo. We find that satellite galaxies lose their gas much more rapidly than their other mass components, including gas located in their central regions. We also find that, although stripping is strongest in the outer regions, the intermediate and inner parts of satellites experience mass loss at early times (i.e. before the galaxy is formally classified as a satellite).

The results of this work quantify and characterise the evolution of both the amount and the spatial distribution of mass in satellite galaxies and complement those of Contreras et al. (2026), which measure the impact of baryons on the positions and velocities of subhaloes. The trends identified in these studies can be incorporated into both empirical galaxy–halo connection models, such as HOD and SHAM-like models, and galaxy formation models built on gravity-only simulations, such as semi-analytical models (SAMs; Kauffmann et al. 1993; Cole et al. 1994; Somerville & Primack 1999; Bower et al. 2006; Lagos et al. 2008; Somerville et al. 2008; Benson 2010, 2012; Jiang et al. 2014; Croton et al. 2016; Lagos et al. 2018; Henriques et al. 2020) and other empirical approaches (Yang et al. 2003; Zheng et al. 2009; Cacciato et al. 2013; Guo et al. 2016; Behroozi et al. 2019; Moster et al. 2020) to improve their description of satellite galaxies. These results also help model the baryon distribution in the context of satellite galaxies, thereby affecting lensing statistics.

The outline of this paper is as follows. Sect. 2 presents the MTNG simulation and describes the galaxy and subhalo samples used in this work. Sect. 3 presents the evolution of the baryonic content of satellite galaxies, and Sect. 4 shows the evolution of their mass profiles. We conclude in Sect. 5. Appendix A presents a parametrisation of the gas content of the central galaxies, which can be combined with the results of this paper to predict the gas content of satellite galaxies. Appendix B examines the intrinsic differences between the central galaxies and satellites shortly after infall, while Appendix C explores additional parametrisations of galaxy properties that may facilitate their modelling in gravity-only simulations. Finally, Appendix D provides an approximate assessment of the bias in satellite gas masses introduced by the subhalo-finding procedure used in MTNG through a comparison with TNG300.

Unless otherwise stated, masses are given in units of h⁻¹ solar masses (h⁻¹ M_⊙) and distances in units of h⁻¹ megaparsecs (h⁻¹ Mpc). All logarithms are in base 10. The virial mass definition used in this work is m_200c, defined as the total mass of a halo enclosed within a sphere whose mean density is 200 times the critical density of the Universe. Throughout this work, the term satellite galaxy refers to the full subhalo, including the associated dark matter subhalo and full gas content, rather than only the central core composed mainly of stellar mass and cold gas.

2. Numerical simulations and galaxy samples

In this section, we present the MTNG hydrodynamic simulation used in this work (Sect. 2.1) and describe the criteria used to select galaxies in each sample for comparison (Sect. 2.2). These samples define the baseline galaxy selections used throughout the analysis.

2.1. The MTNG simulation

To quantify the evolution of satellite galaxy properties, we require a hydrodynamic simulation with a large volume to ensure robust statistics, as well as sufficiently high resolution to resolve multiple satellites per halo even after they have been stripped of most of their mass. One of the few simulations currently available with these characteristics is the largest hydrodynamic run in the MillenniumTNG suite of simulations (Pakmor et al. 2023; Hernández-Aguayo et al. 2023; Barrera et al. 2023; Hadzhiyska et al. 2023b,a; Delgado et al. 2023; Kannan et al. 2023; Bose et al. 2023; Contreras et al. 2023b). This suite comprises tens of gravity-only (usually referred to as dark-matter-only, DMO) and hydrodynamic simulations, together with light-cone galaxy catalogues generated using SAMs. The project extends the 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 over much larger cosmological volumes than previously possible.

The main simulation used in this work is the largest hydrodynamic simulation in this suite, which we refer to as the MTNG simulation. It has a periodic volume of (500 h⁻¹ Mpc)³ and contains 4320³ dark matter particles and an equal number of gas cells, implying an average gas-cell mass of 2.00 × 10⁷ h⁻¹ M_⊙ (and a similar initial mass for newly formed stellar particles) and a dark matter particle mass of 1.12 × 10⁸ h⁻¹ M_⊙. The simulation adopts a cosmology consistent with Planck Collaboration XIII (2016)¹, 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 root-mean-square values (Angulo & Pontzen 2016), a procedure that substantially reduces the impact of cosmic variance, at least for second-order statistics. The simulation was evolved with the moving-mesh code AREPO (Springel 2010), which includes star formation, radiative cooling, the growth of supermassive black holes, and the associated feedback from both supernovae and active galactic nuclei (AGN). The haloes in this simulation were identified using a FRIENDS-OF-FRIENDS (FoF) group finder and a modified version of the SUBFIND-HBT substructure finder (Springel et al. 2021), adapted from the Gadget4 code to AREPO. In this implementation of the substructure finder, gas-cell membership in satellite subhaloes is assigned using previous subhalo membership to find possible candidates, analogously to dark matter and star particles. Because gas cells in AREPO are not strictly Lagrangian and can exchange mass across cell interfaces, this procedure can lead to a systematic underestimation of the diffuse gas mass associated with satellites. This limitation should be taken into account when interpreting the gas content of satellite subhaloes in MTNG.

Galaxies form as agglomerations of star particles whose properties can be measured directly from the simulation. Our simulation includes 265 snapshots, spanning from z = 63 (when the initial conditions were generated) to z = 0. The large number of snapshots provides a unique opportunity to follow in detail the infall of satellite galaxies, during which the relevant physical processes occur on short timescales. However, only the group and subhalo catalogues were available for most snapshots, while particle data were accessible only for a limited number of key outputs. In this work, we used the particle data at z = 0 to study the mass profiles of satellite galaxies.

In Appendix B, we discuss matching the satellite galaxies in the MTNG simulation to the subhaloes of a gravity-only simulation from the MillenniumTNG suite with the same number of dark matter particles, cosmology, initial conditions, and volume, and run with the Gadget-4 code (Springel et al. 2021). We refer to this simulation as MTNG-DMO. To match the catalogues, we used the procedure described in Contreras et al. (2026), which traces satellites back to the snapshot at which they were still centrals using the merger trees of the simulations, and then matches the objects based on their mass and position.

2.2. Target galaxy samples

In this work, we used satellite galaxies selected at z = 0 and traced them back through the merger trees to z ∼ 13. We studied the evolution of four properties: stellar mass, gas mass, r-band magnitude, and U-band magnitude. The synthetic magnitudes used in this work were taken from the MTNG group catalogues and were computed by summing the luminosities of all stellar particles assigned to each subhalo. The luminosity of each stellar particle was computed from stellar population synthesis tables (Bruzual & Charlot 2003) using its age, metallicity, and initial mass, and then convolved with the corresponding photometric filters. The U-band magnitudes are given in the Vega system, while the SDSS r-band magnitudes are given in the AB system. Since our analysis is based on the evolution of each magnitude relative to its value at infall, this choice of absolute photometric zero-point does not affect the trends studied here. We did not apply any additional dust attenuation correction. Galaxies were selected within narrow ranges of host halo mass (m_200c), defined as the mass enclosed within a sphere whose mean density is 200 times the critical density of the Universe, and peak subhalo mass (m_peak), defined as the maximum mass that a subhalo reaches during its evolutionary history. In particular, we selected galaxies in three peak-subhalo-mass bins spanning 10¹¹−10^12.5 h⁻¹ M_⊙, each 0.5 dex wide, and in three host-halo-mass bins spanning 10^12.75−10^14.25 h⁻¹ M_⊙, also 0.5 dex wide. A visual representation of the selected samples is shown in the upper panel of Fig. 1. The coloured lines show the mean number of satellite galaxies per halo as a function of halo mass for the m_peak-selected samples. The vertical shaded regions represent the additional halo-mass selection applied to the galaxy samples. Each sample should therefore be understood as the subset of the m_peak-selected galaxies that also lies within the corresponding halo-mass bin.

Fig. 1. Satellite HOD and stellar-mass relation. Top panel: The solid lines show the satellite HOD for galaxies in the MTNG simulation selected by peak subhalo mass, i.e. the maximum subhalo mass reached by each galaxy during its evolutionary history, while the dashed black line shows all subhaloes with peak subhalo mass above 1011 h−1 M⊙. The vertical coloured bands correspond to the halo-mass selection used throughout this work. Bottom panel: The solid line shows the median stellar mass of satellite galaxies as a function of peak subhalo mass, while the shaded regions represent the 16th and 84th percentiles. The vertical coloured bands correspond to the peak-subhalo-mass selection used throughout this work.

The selection of galaxies as a function of m_peak is intended to avoid biasing the sample by the galaxies’ evolutionary stage. As we show in the next section, galaxy properties such as stellar mass or magnitudes evolve after infall. In addition, over most of the range considered, stellar mass strongly correlates with m_peak, as shown in the bottom panel of Fig. 1, which presents the stellar mass of satellite galaxies in the MTNG simulation as a function of m_peak.

3. The evolution of satellite galaxies

3.1. Evolution of the baryonic content in satellites

We examined the evolution of the stellar mass, gas mass, and r- and U-band magnitudes of satellite galaxies. Figure 2 shows the evolution of these properties for a representative satellite galaxy with peak subhalo mass of m_peak ∼ 10^12.3 h⁻¹ M_⊙ orbiting within a halo of mass m_200c ∼ 10^13.5 h⁻¹ M_⊙. In addition to the four previously mentioned properties, we also show the evolution of the star formation rate and the total subhalo mass of this galaxy, which are often used to characterise satellite evolution. Each dot represents the position of the satellite relative to the halo centre, which is located at the centre of each panel. Each column shows the evolution of a different property, as labelled, with the area of each dot representing the value of the corresponding quantity. The figure shows that properties such as stellar mass and M_r change only modestly over most of the satellite’s lifetime, with stronger evolution occurring at later stages. This is because stellar mass is concentrated in the inner regions of the satellite and is strongly gravitationally bound, which largely protects it from the stripping processes experienced by the subhalo as it moves through the host halo’s hot gas and tidal field. It is expected that properties that are closely linked to stellar mass, such as red-band magnitudes, will also follow a similar evolutionary trend. By contrast, the gas and dark matter are less tightly bound and are more spatially extended, causing them to lose most of their mass during the first pericentric passages, when the tidal field experienced by the satellite is strongest (e.g. Stücker et al. 2023; Aguirre-Santaella et al. 2023). The gas content declines more abruptly than the total mass, with virtually no gas left bound to the satellite before the galaxy completes its second orbit around the host. Finally, we note that the star formation rate also declines rapidly, but exhibits a strong starburst episode around the time of pericentric passage, clearly visible in the Y–Z plane (middle row). Here, ram-pressure compression could temporarily enhance star formation before quenching. This effect is also reflected in the evolution of M_U, which shows a temporary decrease (i.e. an increase in luminosity), since the U band is sensitive to both the stellar mass of the galaxy and its star formation rate.

Fig. 2. Illustration of a satellite with mstell ∼ 1010.5 h−1 M⊙ infalling into a halo of ∼ 1013.5 h−1 M⊙ at z = 0. The solid circles show the satellite’s position relative to the halo centre. The size of each dot scales with the stellar mass, gas mass, star formation rate, and subhalo mass (in the first, second, fifth, and sixth columns, respectively), and with the difference in the r- and U-band magnitudes (in the third and fourth columns), relative to their values at the moment of infall. The larger dotted circle represents one virial radius (r200c).

We aimed to quantify the evolution of larger galaxy samples. To do so, we used the criteria described in Sect. 2.2, which select galaxies according to their host halo mass and peak subhalo mass. Figure 3 shows the evolution of the stellar mass, gas mass, and r- and U-band magnitudes (first, second, third, and fourth rows, respectively) as a function of four quantities: the ratio of the current subhalo mass to its infall value, the pericentric distance normalised by the comoving virial radius of the host halo, the time since infall, and the ratio of the current maximum circular velocity to its infall value. We defined $v_{\max }\equiv \max \sqrt{GM(<r)/r}$ as the total mass in the central regions of subhaloes, which should hence be interpreted as a subhalo property rather than a direct measure of galaxy kinematics. These quantities were chosen because they characterise different aspects of satellite evolution and are commonly available in both hydrodynamic and gravity-only simulations.

Fig. 3. Evolution of the stellar mass, gas mass, and r- and U-band magnitudes of satellite galaxies at z = 0 respectively shown as a function of the remaining subhalo mass fraction (first column), pericentric distance (second column), time since infall (third column), and the remaining maximum circular velocity fraction (fourth column). The different colours represent different host halo masses, while the different line styles represent different peak subhalo masses, as labelled. For simplicity, we show different halo masses only for the intermediate-mpeak sample. The shaded region represents the 16th and 84th percentiles of the distribution.

The remaining subhalo mass fraction can be readily obtained in simulations with merger trees and is expected to correlate with the gas mass of galaxies, since both the gas and the dark matter are spatially extended throughout the subhalo. This correspondence is not exact, however, because gas stripping is additionally affected by ram pressure. This quantity has been used to study the evolution of stellar mass in satellite galaxies (e.g. Smith et al. 2016; He et al. 2026) and to model it in semi-empirical models such as EMERGE (Moster et al. 2020) and in early versions of the empirical model SHAMe (Contreras et al. 2021b).

The pericentric distance, defined as the minimum distance of the galaxy from the halo centre, has been suggested to be related to episodes of rapid mass loss, both in total mass and gas mass, as well as to quenching processes, since tidal forces are expected to be strongest near pericentre (e.g. Gnedin et al. 1999; Wright et al. 2022; Stücker et al. 2023). This picture is consistent with the evolution of the individual galaxy shown in Fig. 2. To facilitate the stacking of many satellites, we normalised this distance by the comoving virial radius of the host halo. The pericentric distance is relatively straightforward to compute in any simulation with merger trees, but its value will always be an upper limit to the true value because of the finite number of snapshots, unless higher-order interpolation between snapshots is performed (for example, in spherical coordinates). This can become a significant limitation in simulations with poor temporal resolution, but not in MTNG, where, as shown in Fig. 2, the large number of snapshots makes this calculation sufficiently accurate.

The time since infall is perhaps the most direct way to characterise the evolution of a satellite galaxy. Although its precision also depends on the temporal resolution, the associated error is easier to estimate than that of the pericentric distance, which is sensitive to the orbital speed of the satellite galaxy. This quantity is also used by the latest version of the empirical model SHAMe to model the fading of satellite galaxies (Contreras et al. 2023a).

The final quantity we used to characterise satellite evolution is the remaining maximum circular velocity fraction (v_max/v_{max, infall}). It is well known that v_max is strongly correlated with stellar mass, and for this reason it is commonly used in the construction of mock catalogues within the SHAM framework (e.g. Chaves-Montero et al. 2016). Unlike the subhalo mass fraction, v_max is primarily sensitive to the central regions of the subhalo and therefore requires higher resolution to be measured accurately, since a larger number of particles per subhalo is needed to achieve good precision. This quantity is used by semi-empirical models such as UniverseMachine (Behroozi et al. 2019) to describe the evolution of satellite galaxies. An additional limitation of using v_max as a proxy for satellite evolution is that, in hydrodynamic simulations, stellar particles eventually contribute to the calculation of v_max, artificially enhancing part of the correlation between these two properties. As previously shown when matching hydrodynamic and gravity-only simulations, a tight correlation between these quantities is expected (Chaves-Montero et al. 2016). Moreover, as shown in Fig. 3, the stellar mass of galaxies does not follow a monotonic relation with v_max, so we do not expect this contamination to dominate the trends presented here. Nevertheless, we warn the reader that this is an additional limitation of using v_max.

The way galaxy properties evolve varies from case to case. The stellar mass of satellite galaxies increases significantly after infall, by around 30%–50%, depending on both the mass of the original galaxy and that of the host halo. By contrast, the gas mass decreases much more rapidly than the stellar mass changes. This behaviour is not captured when the evolution is parametrised in terms of v_max, since a large fraction of the gas is lost while v_max remains nearly constant, thereby washing out the trend. The magnitudes in both bands evolve in all the cases we consider and follow broadly similar trends to each other, although the U-band magnitude fades more rapidly (i.e. M_U increases faster) than the r-band magnitude. We note that the evolution of the magnitudes is not expected to follow the stellar-mass evolution in a one-to-one manner. Although red-band luminosities have a strong correlation with stellar mass in a statistical sense, their evolution also depends on the stellar population age, metallicity, and star-formation history. This dependence is particularly relevant for satellite galaxies, for which environmental processes can suppress star formation and modify the stellar mass-to-light ratio without producing an equivalent change in the total stellar mass. For this reason, we model the stellar-mass and magnitude evolution separately.

We find that all the subhalo properties considered can reasonably capture the evolution of the galaxy properties studied. That said, some appear to be more suitable than others for describing the evolution of specific galaxy properties. In the case of stellar mass, only the change in v_max is able to trace the decrease in stellar mass down to values below 50% of those at infall, while doing so with only weak dependence on galaxy sample and host halo mass. As mentioned above, this quantity strongly correlates with a galaxy’s total stellar mass, as it is determined by the mass distribution in the subhalo’s central regions. When the stellar mass evolution is characterised by this subhalo property, we find only a weak dependence on host halo mass, which facilitates the use of these trends in models of satellite galaxy evolution. A mild dependence on host halo mass is seen only for the lowest-m_peak sample, where the trend deviates for galaxies in the most massive host haloes. In this subhalo sample, galaxies can contain as few as 20 stellar particles (∼ 4 × 10⁸ h⁻¹ M_⊙), meaning that the results are strongly affected by resolution limits. As discussed above, the use of v_max can itself be limited by the resolution of the simulation, which may explain why this disagreement is less evident when the evolution is characterised in terms of the other subhalo properties. The other samples with the same m_peak but different host halo masses show similar results. For this reason, we do not investigate this difference further, but nevertheless show the results to make the reader aware of the limitations of using v_max as a proxy for galaxy evolution. In addition to the maximum circular velocity, changes in the subhalo mass also capture the evolution of stellar mass, although less smoothly than v_max.

For the gas mass, we find that the best way to characterise its evolution is by using the remaining subhalo mass fraction. When expressed as a function of this property, the relation shows only a small scatter and a weak dependence on the galaxy sample chosen. These characteristics make the remaining subhalo mass fraction an ideal candidate for modelling this galaxy property. Small differences between galaxy samples appear only once the galaxy has lost more than 80% of its gas mass. Another subhalo property capable of characterising this evolution is the time since infall, although it shows larger scatter and a stronger dependence on the galaxy sample. We do not find that either pericentric distance or v_max is able to capture the evolution of this galaxy property. In the case of pericentric distance, the evolution of individual galaxies shows that satellites begin to lose gas soon after becoming satellites. Since galaxies become satellites at different virial radii, owing to differences in the shape and mass of the host halo, it is difficult to model a single overall trend for all galaxies. This limitation, however, does not preclude the usefulness of pericentric distance for modelling the evolution of satellite properties but rather indicates that it becomes more difficult to characterise a heterogeneous galaxy sample drawn from different host haloes. In the case of v_max, changes in this quantity occur only after galaxies have already lost a considerable amount of gas and total mass, which makes it unreliable as a proxy for modelling the evolution of this property. We also examined the evolution of the total baryonic content (stellar mass plus gas mass), which evolves similarly to gas mass but with more universal trends (i.e. no significant differences between galaxy samples) when expressed as a function of subhalo mass (not shown here). As discussed in Appendix D, the diffuse gas mass assigned to satellite subhaloes in MTNG may be systematically underestimated. The absolute gas-mass scale should therefore be interpreted with caution. Nevertheless, comparison with TNG300 (Nelson et al. 2018; Springel et al. 2018) suggests that the overall trend is preserved, so this relation still provides a useful description of the relative evolution of satellite gas content.

All the subhalo properties considered can describe the evolution of the magnitudes. When expressed as a function of the remaining subhalo mass fraction, the magnitudes show the greatest scatter among the different parametrisations but still provide a good description of their evolution, with galaxies becoming gradually fainter as the mass fraction decreases. When expressed as a function of pericentric distance, the magnitudes remain roughly constant out to about half the comoving virial radius and then begin to evolve smoothly. When expressed as a function of time since infall, galaxies tend to maintain their magnitudes for 1–2 Gyr and then become fainter, with their magnitudes changing linearly with time. Finally, when expressed in terms of the v_max ratio, galaxies exhibit almost linear evolution with low scatter and show no major differences across different galaxy samples. Although the evolution of M_r and M_U differs in detail, their evolutionary trends are similar, which facilitates their modelling.

We examined the evolution in other photometric bands and find qualitatively similar trends. Bluer bands are more sensitive to recent star formation, as blue, luminous, high-mass stars evolve off the main sequence much more rapidly than their redder, fainter, lower-mass counterparts. As the gas content declines and star formation is suppressed, galaxies become progressively fainter and redder, giving rise to the trends identified in this work. A more detailed discussion is provided by Chaves-Montero & Hearin (2020, 2021).

We also examined the evolution of satellite properties at z = 1 (Fig. 4). We find trends similar to those at z = 0, with the main exception of those expressed as a function of time since infall. In addition, the stellar mass tends to increase slightly more at z = 1 than at z = 0. Despite this difference, the overall agreement between the two redshifts suggests that these trends may be largely universal and could therefore be used in models of satellite evolution independently of the infall redshift of the galaxies.

Fig. 4. Evolution of satellite galaxy properties at z = 1. This figure is similar to Fig. 3, but for z = 1.

Regarding the scatter in these relations, as discussed in Appendix B, the stellar and gas masses of satellite galaxies are correlated with the large-scale environment, which likely contributes to the scatter in these relations. This scatter is also expected to depend on other secondary halo properties, such as concentration or halo age, which may influence the evolution of both central and satellite galaxies (e.g. Zehavi et al. 2018; Artale et al. 2018). A more detailed analysis of the physical origin of this scatter could further improve our understanding of satellite evolution and will be addressed in future work.

3.2. Modelling the evolution of satellite galaxies

Most of the relations presented in the previous section can be used to some extent to model the evolution of satellite galaxies. Although some trends show lower scatter, weaker dependence on the galaxy sample or host halo mass, and are easier to model, no single proxy provides an optimal description of all galaxy properties. For this reason, we present the full set of relevant relations identified in this work, allowing the choice of proxy to depend on the specific scientific goal. To facilitate the modelling of these properties, we provide a parametrisation for each galaxy property that can be used to describe satellite evolution.

These parametrisations can be combined with known relations between galaxy and halo properties for central galaxies, providing an initial value for each galaxy property. Using these initial values, the subsequent evolution can be modelled as satellites fall into their host haloes. For instance, SHAM has been shown to perform well when associating halo properties with galaxy properties, particularly for quantities such as stellar mass or M_r. In the case of M_U, a more complex abundance-matching prescription may be required, such as that proposed by Ortega-Martinez et al. (2024, 2025). Appendix A provides a parametrisation for the gas mass that can be used to assign an initial value to central galaxies by characterising their gas mass as a function of v_max.

The key relations we parametrised are the evolution of gas mass as a function of the remaining subhalo mass fraction, and the evolution of stellar mass and magnitudes as a function of the remaining maximum circular velocity fraction. These relations, together with their corresponding parametrisations, are shown in Fig. 5. For simplicity, we show only the results for the intermediate m_peak sample (11.5 < log(m_peak/h⁻¹ M_⊙) < 12.0). Other galaxy samples show trends similar to those shown here. The decrease in the gas mass was modelled as a shifted and rescaled rational function:

$$\begin{array}{c}\hfill \frac{m_{\mathrm{gas}}}{m_{\mathrm{gas},\mathrm{infall}}}=max(0,\frac{x}{x+A{(1-x)}^{\beta }}(1+\gamma )-\gamma ),\end{array}$$(1)

Fig. 5. Evolution of the stellar mass as a function of vmax/vinfall (top left panel), the gas mass as a function of msub/minfall (top right panel), and the r- and U-band magnitudes as a function of vmax/vinfall (bottom left and right panels) for galaxies with 11.5 < log(mpeak/h−1 M⊙) < 12.0. The symbols represent the data from the MTNG simulation (as shown in Fig. 3), while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled.

where x  ≡  m_sub/m_infall, and A, β, and γ are free parameters. A simpler parametrisation, x/[x + A(1 − x)^β], is also able to reproduce this trend, although with some limitations at gas fractions below 0.2. As a reference, the best-fitting parameters for the log(m_200c/h⁻¹ M_⊙) = 13.5 sample are A = 6.0 ± 0.4, β = 1.17 ± 0.03, and γ = 0.021 ± 0.006.

We emphasise, however, that this parametrisation was calibrated to the gas masses assigned to satellite subhaloes in MTNG. As discussed in Appendix D, its absolute normalisation should be interpreted with caution, although comparison with TNG300 suggests that the overall shape of the relation remains reasonably robust.

For the stellar mass, we modelled its evolution as a power law with an exponential decay:

$$\begin{array}{c}\hfill \frac{m_{\mathrm{stell}}}{m_{\mathrm{stell},\mathrm{infall}}}=\left(\frac{A}{A+1-x}\right)e^{-B(1-x)},\end{array}$$(2)

where x  ≡  v_max/v_infall, and A and B are free parameters. We tested other, more complex parametrisations, but they did not yield significant improvements unless the number of free parameters was increased substantially. As a reference, the best-fitting parameters for the log(m_200c/h⁻¹ M_⊙) = 13.5 sample are A = 16.7 ± 5.1 and B = 1.22 ± 0.19. Other parametrisations available in the literature (e.g. Smith et al. 2016; He et al. 2026) focus primarily on the decay of stellar mass, typically modelled as an exponential decline, and do not capture the increase in stellar mass that we find after infall.

Finally, we modelled the change in magnitude as a linear decay with a negative power law:

$$\begin{array}{c}\hfill \mathrm{\Delta }M=\mathrm{M}-{\mathrm{M}}_{\mathrm{infall}}=-A{(1-x)}^B,\end{array}$$(3)

where x  ≡  v_max/v_infall, and A and B are free parameters. We also tested an even simpler form, a linear decay ΔM = −A(1 − x), which provides a reasonably good description over the full v_max range. As a reference, the best-fitting parameters for the log(m_200c/h⁻¹ M_⊙) = 13.5 sample are A = 4.8 ± 0.1 and B = 0.85 ± 0.02 for the r band, and A = 7.8 ± 0.8 and B = 0.71 ± 0.06 for the U band. All fits assume a constant uncertainty for each data point, scaled such that the best fit yields a reduced chi-squared equal to one. The best-fitting parameters for all the main z = 0 samples, together with the scatter and RMS residuals of the fits, are provided in Appendix E.

The choice of the best subhalo parameter for modelling these galaxy properties was based on the overall performance and ease of use described in the previous section. Nevertheless, the most suitable subhalo property for characterising a given galaxy property may depend on the specific scientific application. In Appendix C, we examine the performance of the parametrisations described in this section when applied to other subhalo properties, finding that they can reproduce the evolution of galaxy properties in most cases. We also provide an additional parametrisation valid across all halo masses, and show that these relations remain valid at z = 1. These results should facilitate the implementation of analytic models to follow satellite evolution in gravity-only simulations across a broader range of scientific applications. The additional parametrisation also provides a better fit to the stellar-mass evolution, reducing some of the residuals seen in Fig. 5. In particular, the simple model does not reproduce all the details of the relation, and small systematic residuals remain at high and low values of v_max/v_infall. We nevertheless kept the simpler parametrisation as our default model, since the more flexible form substantially increases the number of free parameters, is not forced to recover the correct normalisation at infall, and is more likely to capture simulation-specific features. The simple model should therefore be interpreted as a compact description of the median trend as opposed to a high-precision fit to all the features of the MTNG relation. A more detailed analysis of the dependence of these trends on peak subhalo mass and host halo mass, and of the dependence of higher-order parameters on the scatter of these relations, could provide even more accurate predictions; however, this lies beyond the scope of the present paper.

4. The evolution of the mass profile of satellite galaxies

In the previous section, we quantified and modelled the evolution of satellite galaxy masses and magnitudes as they fall into their host haloes. We now turn to the evolution of the mass profiles of their different components. Since the binding energy is highest in the central regions of galaxies, mass loss is expected to begin in the outer parts of their subhaloes, primarily removing dark matter and gas, and only affecting the stellar component once the galaxies have already lost most of their mass (e.g. Cortese et al. 2021). This picture is consistent with the results in Sect. 3.1, with the exception that galaxies lose most, if not all, of their gas much faster than they lose their dark matter. We expect the diffuse outer gas to be stripped more efficiently than the dark matter because, although both components are subject to tidal stripping (Tormen et al. 1998; van den Bosch et al. 2018), only the gas is affected by ram pressure, which can efficiently remove weakly bound material once the external pressure exceeds the local gravitational restoring force (McCarthy et al. 2008, see also Ayromlou et al. 2021).

In Fig. 6, we show the evolution of the mass profiles of satellites with 11.5 < log(m_peak/h⁻¹ M_⊙) < 12.0 and 12.0 < log(m_peak/h⁻¹ M_⊙) < 12.5 as a function of distance from the galaxy centre, for the dark matter, gas, and stellar mass components at z = 0. We express the distance from the galaxy centre in units of the comoving virial radius of the satellites at the moment of infall, which facilitates stacking galaxies that fell in at different times. Stacking was performed independently of host halo mass, since we did not find significant differences in the mass profiles across halo masses. Although more massive haloes provide a more hostile environment, selecting satellites by their m_sub/m_peak ratio allows us to compare galaxies at a similar evolutionary stage, largely independently of host halo mass. We do not include the lowest-m_peak galaxy sample because the number of particles per object is low, especially for the stellar mass component. The y-axis represents the total mass in a spherical shell (i.e. we do not normalise by the shell volume). The different coloured lines represent the m_sub/m_peak bins. This ratio is similar to m_sub/m_infall and allows us to characterise the evolution of the profiles independently of the definitions of central and satellite galaxies. We group galaxies in bins of 0.2 in m_sub/m_peak (as labelled). We also include a galaxy sample with m_sub/m_peak > 0.97, which represents satellites shortly after infall, as well as the mass profile of central galaxies.

Fig. 6. Mass profiles of satellite subhaloes with peak subhalo masses between 1011.5 and 1012 h−1 M⊙ (top row) and between 1012 and 1012.5 h−1 M⊙ (bottom row). The first, second, and third columns show the subhalo mass, gas mass, and stellar mass profiles, respectively. The solid coloured lines represent different levels of subhalo mass loss, while the dashed black lines represent the profiles of central subhaloes. Each profile is normalised by the comoving virial radius (r200c) of the satellite at the moment of infall.

In agreement with the results in Sect. 3.1, we find that the gas component is the first to be removed, initially from the outer regions and later from the inner ones. The dark matter component shows a similar behaviour, but on a slower timescale. The stellar mass component first increases and only later begins to decrease, and it is not yet clear whether this stripping proceeds from the outside in or follows a more complex behaviour. These trends are shown more clearly in Fig. 7, where we present the ratio of the different subhalo-mass-ratio samples to the 0.8 < m_sub/m_peak < 1 sample. We do not use the sample of satellites shortly after infall for the denominator, since it is noisier because of its limited number of galaxies. It is important to note that, although most of the mass is lost from the outer regions, there is a continuous loss of gas in the central regions. This suggests that the simple picture in which mass is removed smoothly from the outside in provides a useful approximation in many cases but does not fully capture the complexity of satellite evolution. The figure may also indicate that dark matter is stripped from the central regions, although this could partly reflect the fact that the dark matter mass profiles do not fully converge to unity on the scales shown here (see Stücker et al. 2023 for a discussion of this issue).

Fig. 7. Ratios of satellite mass profiles. This figure is similar to Fig. 6, but shows the ratio of the profiles to the 0.8 < msub/mpeak < 1 sample.

When comparing the evolution of the two satellite samples selected by m_peak, we find overall qualitative agreement in how galaxies lose dark matter and gas but not in how their stellar mass component evolves. Lower-mass galaxies show a more pronounced increase in stellar mass, peaking at around 10% of their virial radius. Only when the subhalo mass has fallen significantly does the stellar mass begin to decrease. The more massive sample loses mass first from the outer regions and later from the inner ones, while the less massive sample shows no clear trend. The overall good agreement between the results for the two galaxy samples, except for the stellar mass component, which contributes only a small fraction of the total subhalo mass, suggests that the way satellite galaxies lose mass may be largely independent of the masses of the satellite and the host halo. Unfortunately, given the limited number of samples that can be studied because of the resolution and volume limitations and the fact that we can examine only one galaxy formation prescription, we cannot generalise this result to a universal relation. A more detailed study of satellite evolution in other simulations is left for future work.

In addition to the evolution of satellite galaxies, Fig. 6 also reveals differences between the profiles of satellites shortly after infall and central galaxies. These results are independent of the normalisation adopted for the distance to the galaxy centre, as well as of the bin width in m_sub/m_peak. These differences arise because central galaxies are intrinsically different from satellites, even shortly after infall, with centrals having more gas and less stellar mass, consistent with the results of Contreras et al. (2015) for SAMs and He et al. (2026) for hydrodynamic simulations. This may be partly related to the fact that, at z = 0, more concentrated haloes tend to reside in denser environments, an effect known as halo assembly bias (e.g. Gao et al. 2005; Gao & White 2007), and these environments are also known to host more massive and less star-forming galaxies (e.g. Zehavi et al. 2018; Artale et al. 2018). By selecting massive satellite galaxies, as done in this work, we selected objects in denser environments than those associated with a less-biased sample of central galaxies. The differences between centrals and satellites are particularly pronounced for the gas component. Part of this difference also arises because satellites reach their peak gas mass 0.5–1 Gyr before they reach their peak subhalo mass, and up to 2 Gyr before they become satellites. This suggests that the infall process, and therefore gas stripping, begins earlier than implied by the moment at which the galaxy is formally classified as a satellite, while the subhalo may still increase in mass, probably owing to interactions with the outer regions of the host halo it will eventually merge into. Since gas loss is abrupt, with galaxies losing most of their gas within 2–3 Gyr after infall, satellites may already have lost a considerable fraction of their gas before they are classified as satellites or before they reach their peak subhalo mass, which may explain the behaviour seen in Figs. 6 and 7. This result suggests that, for some applications, the usual central–satellite classification may not be ideal. Empirical models such as HODs rely entirely on this classification and generally do not distinguish between the overall central population and centrals that are about to become satellites. Our results therefore suggest that a more conservative definition of satellite galaxies, in which systems are classified as satellites at earlier times, may be required. Although the focus of this paper is not the intrinsic differences between central and satellite galaxies, Appendix B includes a more detailed analysis of the origin of the differences discussed in this section.

5. Summary

In this work, we examined the evolution of the masses and mass profiles of satellite galaxies in the cosmological hydrodynamic MTNG simulation. We first identified satellite galaxies at z = 0 and traced them back to z ∼ 13. We quantified the evolution of their dark matter mass, stellar mass, gas mass, and their r- and U-band magnitudes as a function of the ratio of the current subhalo mass to its value at infall (m_sub/m_infall), the pericentric distance, the time since infall, and the ratio of the current maximum circular velocity to its value at infall (v_max/v_infall) (Fig. 3). All of these quantities can be obtained from gravity-only simulations with merger trees and can therefore be used to model the evolution of galaxy properties in these simpler simulations. The main findings of this part of the work are as follows:

• The stellar mass increases as the satellite evolves, and only begins to decrease once the galaxy has lost a substantial fraction of its subhalo mass or maximum circular velocity. We find that the remaining subhalo mass fraction and the maximum circular velocity ratio can both be used to characterise stellar mass evolution.

• The gas mass decreases abruptly as the satellite evolves, with most of the gas mass lost once the galaxy has lost only about half of its total subhalo mass. The remaining subhalo mass fraction and the time since infall provide the best characterisation of this evolution, with the former showing lower scatter. As discussed in Appendix D, the absolute gas-mass scale should nevertheless be treated with caution, as the gas masses assigned to satellite subhaloes in MTNG may be biased low.

• The magnitudes in both bands are reasonably well described by all the tested proxies. Both bands follow similar trends, although the U-band magnitude fades faster than the r-band magnitude.

• The trends found at z = 0 are largely consistent with those found at z = 1, suggesting that they may be broadly universal and therefore applicable to models of satellite evolution across different infall redshifts (Fig. 4).

Most trends show little dependence on the peak subhalo mass of the galaxies or on the mass of their host haloes. The v_max ratio is found to be particularly insensitive to these selections in most cases, with the exception of the gas, for which the remaining subhalo mass fraction appears to be the least sensitive proxy. These results suggest that simple parametric models can capture the main median trends of satellite evolution.

We provided simple parametrisations for the evolution of the baryonic properties studied in this work, with the aim of modelling satellites using empirical and semi-empirical models as well as SAMs. We characterised the evolution of stellar mass and magnitudes as a function of v_max/v_infall, and the evolution of gas mass as a function of m_sub/m_infall. We find that these relations can be described reasonably well with only two or three free parameters, although systematic residuals remain in some regimes (Fig. 5). These parametrisations can be combined with a simple subhalo abundance-matching prescription to assign initial galaxy properties before satellite infall, thus enabling the assignment of galaxy properties to subhaloes in gravity-only simulations. To complement this framework, Appendix A provides a parametrisation for assigning gas mass to central subhaloes in gravity-only simulations, since this property is not well characterised by the standard subhalo abundance-matching formalism. We note, however, that due to the intrinsic differences between standard centrals and centrals shortly before becoming satellites, or satellites shortly after infall (Appendix B), this step may require additional adjustments to perform optimally.

We then focused on the evolution of the dark matter, stellar mass, and gas mass profiles of satellites at different evolutionary stages. We selected galaxies at z = 0 with different ratios of the current subhalo mass to the maximum subhalo mass reached during their evolutionary history (m_sub/m_peak), for two galaxy samples defined by their m_peak values (Figs. 6 and 7). Since we found similar trends across different host halo masses, we stacked all galaxies together to maximise the statistical power of the results. The main findings of this part of the work are as follows:

• In agreement with the results shown in Fig. 3, we find that the gas mass is lost much faster than the total subhalo mass, whereas the stellar mass increases as the satellite evolves and only later begins to decrease, once most of its subhalo mass has been lost.

• The loss of all mass components begins preferentially in the outer regions of the galaxy, and as the satellite evolves, mass loss extends progressively towards the inner regions.

• A significant fraction of the gas mass does not follow this simple trend, with a continuous loss of mass occurring at all radii.

• While the mass profiles differ between the two m_peak samples, the ratio of the profiles of satellites in a late evolutionary stage to those of satellites shortly after infall is similar for the gas and dark matter components, which together account for most of the subhalo mass.

• The stellar-mass profiles show clearer differences between the two m_peak samples, mostly owing to the different starburst episodes experienced by satellites after infall.

The overall agreement in the evolution of the gas and dark matter components suggests, once again, that the internal mass profiles of satellite galaxies could be parametrised relatively easily if needed. That said, the limited number of galaxy samples studied in this work, together with the use of a single galaxy formation model, strongly limits our ability to make stronger claims about the universality of satellite evolution. Nevertheless, the trends identified here suggest that these relations could be recalibrated for different galaxy formation models and used to improve the modelling of satellite galaxies.

In Contreras et al. (2026), we studied how baryons alter the orbits of satellite galaxies, then quantified and parametrised the changes in the positions and velocities of galaxies relative to the subhaloes of gravity-only simulations. The results of the present work extend those findings by providing a more complete picture of how galaxies evolve after becoming satellites. In future work, we plan to incorporate the results of both studies into the extended subhalo abundance matching (SHAMe; Contreras et al. 2021a,b) model to improve clustering predictions in the non-linear regime. In addition, the differences in the mass profiles will be used in the BACCO baryonification model (Aricò et al. 2020, 2021; Aricò & Angulo 2024; Burger et al. 2026) to further improve the precision of baryonification emulators. These results will also be used to extend the range of validity of galaxy–galaxy lensing modelling to smaller scales, such as the GalaxyEmu emulator (Contreras et al. 2023b; Mahony et al. 2026) and the model by Zennaro et al. (2025), as well as to improve the prescriptions for the satellite contribution to the kinetic Sunyaev Zel’dovich signal (see e.g. McCarthy et al. 2025).

Acknowledgments

We thank Jonas Chaves-Montero and Jens Stücker for some useful suggestions and comments. SC acknowledges the support of the “Ramón y Cajal” fellowship (RYC2023-043783-I). SC also acknowledges the support of the “Ayudas para Atracción de Investigadores con Alto Potencial” (2025/00000640) from Universidad de Sevilla. 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]. Author Contributions Statement: The idea for this project was developed by SC, with essential contributions from REA and GA. SC analysed the data, performed the calculations, and wrote the manuscript. REA, GA and LO provided critical comments throughout the development of the project. SB, LH, RP, and VS developed the MTNG simulation and provided feedback on the manuscript.

References

Appendix A: Modelling the gas content of central galaxies

In Sect. 3.2, we present a parametrisation for the evolution of satellite properties as a function of subhalo properties available in gravity-only simulations. These parametrisations can be used by empirical, semi-empirical, and SAMs to improve their description of satellite galaxies. These relations are expressed as the ratios of stellar mass and gas mass to their values at infall, and as changes in the magnitudes. Modelling these evolutions therefore requires an initial value for these properties at the time when the galaxies were still centrals. For stellar mass and the magnitudes, a subhalo abundance-matching formalism can be used, in which a halo property such as halo mass or v_max is rank-ordered and matched to the expected stellar mass or magnitude distribution. Subhalo properties such as peak subhalo mass and peak maximum circular velocity can also be used, although for central galaxies they may not yield a significant improvement in performance. In the case of the U-band magnitude, a more complex abundance-matching prescription, such as that proposed by Ortega-Martinez et al. (2024, 2025), may be required to model this property accurately.

The only property that cannot be fully modelled with a standard abundance-matching formalism is the total gas mass. We expect the total gas mass, normalised by the subhalo mass, to be close to Ω_b/Ω_M, which corresponds to the maximum baryonic mass available within the halo. The difference between this approximation and the true value arises from the stellar-mass and black-hole-mass components, which are expected to be subdominant, and from gas that has been expelled or redistributed by feedback processes. We express this ratio as a function of v_max rather than subhalo mass, since we find that the corresponding relation shows slightly less scatter, although in practice either quantity could be used. We note that related cold-gas scaling relations for central galaxies have been studied in the context of disc gravitational instability (Romeo 2020); Romeo:2020b.

To facilitate the implementation of this relation in galaxy population models, which are usually built on gravity-only simulations, we characterise the gas-to-subhalo-mass ratio as a function of the v_max of the matched haloes in the MTNG-DMO gravity-only simulation. We perform the halo matching between the two simulations by identifying pairs of objects with similar positions and masses. To ensure a clean selection, we use only haloes that are matched bidirectionally, i.e. haloes for which an object in the hydrodynamic simulation is matched to an object in the gravity-only simulation, and for which that same object in the gravity-only simulation is in turn matched back to the original halo in the hydrodynamic simulation.

The gas fraction of central galaxies in the MTNG simulation at z = 0 is shown in Fig. A.1. This relation shows a local peak at about v_max ∼ 200 km s⁻¹, which corresponds to a halo mass of ∼10¹²h⁻¹M_⊙, the scale at which AGN feedback becomes efficient. The gas fraction then decreases before rising again, reaching values close to the cosmological baryon fraction. This shape is similar to that found in other hydrodynamic simulations (e.g. Schaye et al. 2023), suggesting that, although the relation may differ across galaxy formation models, the same fitting function could still be used after recalibrating its free parameters. We parametrise this relation as the sum of a Gaussian, which reproduces the first peak, and an error function, which models the subsequent increase:

$$\begin{array}{c}\hfill f(x)=c+Ag(x)+Bs(x),\end{array}$$(A.1)

Fig. A.1. Gas fraction of central subhaloes in the MTNG simulation, normalised by Ωb/ΩM, as a function of the peak value of the maximum circular velocity reached by each galaxy during its evolutionary history (i.e. vpeak). The solid coloured line represents the median of the distribution, while the shaded region represents the 68th percentile. The dashed black line represents the fit to this distribution (see Appendix A for more details).

where x ≡ log(v_max/km s⁻¹), c is a constant, and g(x) is a Gaussian:

$$\begin{array}{c}\hfill g(x)=exp[-\frac{1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2],\end{array}$$(A.2)

while

$$\begin{array}{c}\hfill s(x)=\frac{1}{2}[1+\text{erf}\left(\frac{x-x_t}{w}\right)]\end{array}$$(A.3)

is the error function. This form successfully reproduces the relation between gas fraction and v_max over the full v_max range. The best-fitting parameters are c, A, μ, σ, B, x_t, and w = 0.167, 0.249, 2.254, 0.122, 0.505, 2.677, and 0.168.

As discussed in Sect. 4 and as discussed in Appendix B, there are intrinsic differences between central galaxies and satellites shortly after infall, with the latter population potentially providing a more appropriate starting point for modelling satellite evolution. We therefore also examine the gas-mass ratio for these galaxies and fit Equation (A.1), finding the best-fitting parameters to be c, A, μ, σ, B, x_t, and w = 0.140, 0.319, 2.248, 0.122, 0.797, 2.68, and 0.132. To match satellites shortly after infall to the dark matter haloes in the gravity-only simulation, we use the method developed in Contreras et al. (2026), which consists of tracing satellites back to the point at which they were still centrals and had a subhalo mass equal to 80% of the peak subhalo mass they would eventually reach, and performing the matching at that redshift.

Appendix B: The differences between central galaxies and satellites shortly after infall

In this appendix, we investigate the origin of the differences in stellar mass and gas mass between standard central galaxies and satellites shortly after infall. In Sect. 4, we find that central galaxies have less stellar mass and more gas than satellites shortly after infall, defined here as satellites with m_sub/m_peak > 0.97. These differences are particularly pronounced for the gas component, for which central galaxies contain substantial amounts of gas in the outer regions of their haloes, whereas satellites shortly after infall lack it.

We quantify these differences in Fig. B.1, where we show the median stellar mass and gas fraction as a function of m_peak. The black shaded region corresponds to the 16th and 84th percentiles of the distribution for central galaxies, and the coloured vertical shaded region corresponds to the m_peak selection used in this work. In terms of stellar mass, satellites shortly after infall have 40%–60% more stellar mass than the average central galaxy; at the same time, they can be up to 60% less massive than the average satellite galaxy at z = 0. For the gas fraction, the relative difference between central galaxies and satellites shortly after infall is about 20%–50%.

Fig. B.1. Stellar and gas masses of centrals and satellites shortly after infall. Top panel: The solid black line and the shaded regions correspond to the median of the distribution for central galaxies and the 16th and 84th percentiles, respectively. The orange line shows the median stellar mass for all satellite galaxies, while the blue line shows only satellites with msub/mpeak > 0.97 (representing satellites shortly after infall). Bottom panel: The same comparison is shown for gas mass normalised by peak subhalo mass. The lower panel in each figure shows the ratio of the coloured lines to the black line (central galaxies).

We identify two main causes for the differences in the properties of these two galaxy populations. The first is the effect of environment: satellite galaxies typically reside in denser environments than central galaxies. The second is that centrals shortly before infall already exhibit an evolutionary behaviour similar to that of satellite galaxies several hundred Myr before they are formally classified as satellites.

B.1. Large-scale environment

Since more massive haloes exhibit a stronger cosmological bias, i.e. they reside in denser environments, and mergers occur more frequently in such environments, it is expected that satellite galaxies also inhabit denser environments than central galaxies of similar mass even before infall. The galaxy population differs across environments even at fixed halo mass, an effect known as occupancy variation. In Zehavi et al. (2018) and Artale et al. (2018), the authors show that galaxies hosted by haloes of the same mass have higher stellar masses in denser environments, consistent with the results found in this paper. This effect is caused by the correlation between large-scale environment and secondary halo properties, such as halo concentration, which are themselves correlated with the galaxies they host; this is commonly referred to as halo assembly bias (Gao et al. 2005; Gao & White 2007); Contreras:2019.

In the context of halo assembly bias, low-peak-height haloes, i.e. low-mass haloes at fixed redshift and cosmology, are expected to be more biased when they are more concentrated. This is shown in Fig. B.2, which presents the halo bias measured in the MTNG-DMO simulation for all haloes, as well as for the upper and lower 20% in concentration, in bins of halo mass. Following Contreras et al. (2021c), we compute the halo assembly bias as the ratio of the halo–matter cross-correlation function for haloes in a given mass bin to the theoretical matter correlation function. To minimise the impact of non-linearities and cosmic variance on the cross-correlation measurement, we compute the clustering on scales between 6 and 20 h⁻¹Mpc and average this ratio over those scales. The cross-correlation function was computed using the public code CORRFUNC (Sinha & Garrison 2019), while the theoretical matter correlation function was obtained from the HALOFIT fitting function (Takahashi et al. 2012), applied to the linear power spectrum from CLASS (Lesgourgues & Tram 2011), and then Fourier-transformed into configuration space.

Fig. B.2. Halo bias in the MTNG-DMO simulation. The halo bias is computed as the ratio of the halo–matter cross-correlation function to the theoretical matter correlation function. The dotted and dashed lines show the biases of the haloes in the lower and upper 20% of the concentration distribution, respectively, in bins of halo mass (i.e. the concentration halo assembly bias). The shaded regions correspond to the bins in mpeak used in this work, which should approximately map onto similar halo-mass ranges.

As in the previous figures, we highlight in Fig. B.2 the halo-mass ranges that approximately correspond to the peak-subhalo-mass cuts used to classify galaxies in this work. While the halo mass in the gravity-only simulation does not map perfectly onto the peak subhalo mass of its hydrodynamic counterpart, we expect only minor differences in the halo selection for central galaxies. A halo with higher concentration is expected to host a central galaxy with higher stellar mass and lower star formation rate, which is related to its gas mass (see Zehavi et al. 2018 for a more detailed discussion). We also note that, in agreement with the differences between central galaxies and satellites shortly after infall found in this work, the differences between the high- and low-concentration samples decrease with increasing halo mass.

To quantify the impact of the large-scale environment on the properties of central galaxies, we show in Fig. B.3 the median stellar mass and gas mass normalised by peak subhalo mass as a function of peak subhalo mass for central galaxies living in the 20% densest and least dense environments. To compute the environment, we use the individual bias-per-object methodology (Paranjape et al. 2018), which effectively measures the bias of an individual object. In agreement with expectations, we find that centrals in denser environments have more stellar mass and less gas than centrals in less dense environments, although the environmental effect is weaker than the difference shown in Fig. B.1. The same conclusion can be reached from Fig. B.4, which shows the mass profiles of centrals living in the densest and least dense environments. Again, centrals in denser environments have more stellar mass and less gas than those in less dense environments. The differences are concentrated mainly around the peaks of each profile and in the outer regions, similarly to the differences found between central galaxies and satellites shortly after infall, although with a smaller amplitude (Fig. 6). It is possible that the environment is not itself the direct cause, but rather traces another effect that drives part of this difference, or that it is only partially responsible for the trends found in this work. That said, among the two baryonic quantities, stellar mass appears to be more strongly affected by environment than gas mass, whereas gas shows larger differences between central galaxies and satellites shortly after infall. This may indicate that environmental effects are not the only drivers of these differences.

Fig. B.3. Stellar and gas masses as a function of environment. Top panel: The solid black line and the shaded regions correspond to the median of the distribution for central galaxies and the 16th and 84th percentiles, respectively. The dashed and dotted lines represent central galaxies in haloes residing in the 20% most and least dense environments, respectively. The environmental classification is based on the individual bias-per-object and is performed in logarithmic bins of halo mass. Bottom panel: The same comparison is shown for gas mass normalised by peak subhalo mass.

Fig. B.4. Mass profiles of centrals in dense and underdense environments. This figure is similar to Fig. 6, but for centrals in the 20% densest and least dense environments. The environmental classification is based on the individual bias-per-object and is performed in logarithmic bins of halo mass.

B.2. Pre-evolution of central galaxies

Another effect responsible for the differences between central galaxies and satellites shortly after infall is the evolution of central galaxies before they become satellites. If centrals begin to be affected by a nearby halo before the time they are formally identified as satellites, it is possible that both their mass components and their internal mass distribution have already changed significantly by the time of infall. To test this hypothesis, we measure the times at which satellite galaxies reach their peak subhalo mass (sometimes referred to as the time of m_peak), stellar mass, and gas mass, and compare them with the time of infall (Table B.1). We find that, for the two most massive samples considered in this work, both the peak subhalo mass and the peak gas mass are reached before the galaxy becomes a satellite. This implies that mass loss begins before the galaxy is formally classified as a satellite. It is important to note that the time at which a galaxy reaches its peak gas mass precedes the time at which it reaches its peak subhalo mass. This picture suggests that, as a galaxy approaches a more massive halo and begins to interact with the outer regions of the host halo, it can still increase its total mass through accretion while its gas is already being stripped.

To build a more general picture of the evolution of centrals before they become satellites, we measure the change in mass as a function of time for both central and satellite galaxies (Table B.2). The mass changes are computed over two six-snapshot windows (∼0.5 Gyr each): one between the time of infall and ∼0.5 Gyr before infall, and another between ∼0.5 and ∼1 Gyr before infall. We find that satellite galaxies show a larger increase in total mass than central galaxies, by about 60%–100%, in both time windows and across both galaxy samples considered. The stellar mass also shows a larger increase, even more pronounced than for the total mass, by about 250%–280%. This suggests that interaction with the nearby halo may trigger the conversion of cold gas into stars, which could also help to explain the differences found between central galaxies and satellites shortly after infall. Finally, for the gas component, we find that its evolution changes significantly as infall approaches. During the period between 1 Gyr and 0.5 Gyr before infall, the increase in gas mass is relatively similar for centrals and satellites, although this still corresponds to only a modest increase for satellites, since during the same period their subhalo masses grow faster than those of central galaxies. In the last 0.5 Gyr before infall, the gas mass barely increases for the more massive satellite sample, while the less massive sample already begins to lose gas. These trends are the same as those found for satellite galaxies elsewhere in this work, strongly suggesting that centrals begin to behave similar to satellites well before they are formally classified as such. We would like to highlight that these mass variations exhibit large scatter and that these numbers should be used only to understand satellite evolution, not for modelling unless a more dedicated study is performed.

We find that both the environmental differences and the distinct pre-infall evolution of centrals shortly before becoming satellites are linked to the differences between central galaxies and satellites shortly after infall. The results presented here, however, cannot disentangle the correlation between these two effects, nor determine whether other processes also contribute to these differences. A more extensive study of this issue would be of clear interest, but lies beyond the scope of this paper. In future work, we plan to investigate the main cause of these differences in more detail. The conclusions of such a study would not only improve our understanding of galaxy formation physics, but also help refine the construction of mock catalogues and potentially clarify some aspects of halo assembly bias.

Appendix C: Additional parametrisations

In this section, we show how the parametrisations presented in Sect. 3.2 perform when other subhalo quantities are used. In Fig. C.1, we show the evolution of the stellar mass and the r- and U-band magnitudes as a function of m_sub/m_infall for galaxies with 11.5 < log(m_peak/h⁻¹M_⊙) < 12.0. We do not include the evolution of the gas mass, since that property is already shown in Fig. 5. The remaining subhalo mass fraction is more commonly used to describe satellite evolution, as it tends to be more readily available than v_max/v_infall. We find that, by using the same parametrisation as for v_max/v_infall, we can recover the evolution of most properties almost perfectly. The only property that presents some difficulty in the modelling is the stellar mass. In this case, we are not able to reproduce the abrupt decrease in stellar mass that occurs when m_sub/m_infall ∼ 0.1. Although such a feature could be expected (Errani et al. 2022), we warn the reader about possible contamination in the calculation of the remaining subhalo mass due to the presence of stellar particles in the core of the galaxy, which are absent in gravity-only simulations.

Fig. C.1. Evolution of stellar mass and magnitudes as a function of remaining subhalo mass fraction. This figure is similar to Fig. 5, but shows the evolution of the stellar mass and the r- and U-band magnitudes as a function of msub/minfall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0. The symbols represent the data from the MTNG simulation, as shown in Fig. 3, while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled.

We then test the parametrisations when modelling the evolution of galaxy properties as a function of time since infall. This quantity does not share the limitation of being computed from the mass distribution within the subhalo, as is the case for m_sub/m_infall and v_max/v_infall, which may differ between hydrodynamic and gravity-only simulations because, in the former, baryonic particles and cells, especially stellar particles, contribute to the inner mass distribution of haloes, whereas this component is absent in gravity-only simulations. We show the evolution of the galaxy properties as a function of time since infall, together with their respective parametric forms, in Fig. C.2. We again use the parametrisations described in Sect. 3.2, although we redefine the fitting variable as x ≡ 1 − t_{since infall}/(10 Gyr). Similar to Fig. C.2, we can successfully model the evolution of the r- and U-band magnitudes, as well as the gas mass. We only struggle to reproduce the final stage of the stellar-mass evolution, for which more complex modelling may be required.

Fig. C.2. Evolution of satellite properties as a function of time since infall. This figure is similar to Fig. 5, but shows the evolution of the stellar mass, gas mass, and r- and U-band magnitudes as a function of time since infall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0. The symbols represent the data from the MTNG simulation, as shown in Fig. 3, while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled.

In this paper, we focus mostly on the evolution of satellite galaxies at z = 0. Nonetheless, Fig. 4 shows that the evolution of galaxy properties at z = 1 is similar to that at z = 0. In Fig. C.3, we show that the parametrisations used throughout this paper remain valid at z = 1 by fitting the evolution of galaxy properties at this redshift. As a reference, we also show the parametrisation of the evolution of galaxy properties for the intermediate-host-halo-mass sample at z = 0. Although the evolution at the two redshifts is not identical, both follow the same trends and remain sufficiently similar to suggest that the parametrisations presented in this work could be used for satellites identified at different redshifts.

Fig. C.3. Analytic fits to satellite galaxy evolution at z = 1. This figure is similar to Fig. 5, but shows the evolution of the stellar mass as a function of vmax/vinfall, the gas mass as a function of msub/minfall, and the r- and U-band magnitudes as a function of vmax/vinfall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0 at z = 1. The symbols represent the data from the MTNG simulation, as shown in Fig. 3, while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled. As a reference, the dotted orange line shows the evolution of the intermediate-halo mass sample at z = 0.

The final set of parametrisations we consider uses the same combination of galaxy and subhalo properties as Fig. 5, but adopts a common parametric form for all host halo masses. For this fit, we add one free parameter to each parametrisation, modifying the values of γ, A, and B in equations (1), (2), and (3), respectively. The new model assumes a linear dependence of these parameters on the logarithm of the halo mass, such that P = P₀ + dP [13 − log(m_200c/h⁻¹M_⊙)], where P₀ is the original free parameter and dP is a new free parameter that captures the dependence on halo mass. As shown in Fig. C.4, all galaxy properties can be well reproduced at all halo masses with this new parametrisation. Although promising, a wider halo-mass range must be explored to demonstrate properly that this parametrisation can be applied across the full halo population of the simulation.

Fig. C.4. Evolution of the stellar mass as a function of vmax/vinfall, the gas mass as a function of msub/minfall, and the r- and U-band magnitudes as a function of vmax/vinfall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0. The different colours represent different host halo masses, as labelled. The symbols represent the data from the MTNG simulation, while the dash-dotted lines represent an analytic model fitted simultaneously to all halo masses.

Regarding the fit to the stellar-mass evolution, we also test symbolic regression in an attempt to provide a more accurate parametrisation of the data, even if it lacks the physical motivation of the form used throughout this work. We find that a parametrisation of the form:

$$\begin{array}{c}\hfill f(x)=\frac{Ax}{{[-exp(Bx)+ln(log(M))]}^{C-Dx}},\end{array}$$(C.1)

where x ≡ v_max/v_infall and A, B, C, and D are free parameters. This parametrisation provides a very close fit to the data, as shown in Fig. C.5, and is particularly useful for reproducing the stellar-mass evolution, which is the most difficult relation to reproduce with the fiducial model. However, its highly flexible form suggests that it may struggle to generalise to other implementations. In addition, because we do not force it to take a value of unity when v_max/v_infall = 1, it may present difficulties when implemented in gravity-only simulations. We nevertheless include this parametrisation to illustrate the potential of symbolic regression to provide highly accurate fits in cases where reproducing a specific model is the main goal.

Fig. C.5. Evolution of stellar mass as a function of vmax/vinfall. The dash-dotted line represents an analytic model, obtained through symbolic regression, fitted simultaneously to all halo masses.

Appendix D: Assessing the impact of the SUBFIND-HBT gas-assignment bias on satellite gas masses

In this appendix, we provide an approximate assessment of the impact of the gas-cell assignment procedure used by the version of SUBFIND-HBT adopted in MTNG. As discussed in Sect. 2, this implementation can systematically underestimate the diffuse gas mass associated with satellite subhaloes. To estimate the magnitude of this effect, we compare the gas-content trends measured in MTNG with those measured in TNG300, which is not affected by this issue in the same way.

The TNG300 simulation has a periodic box of 205 h⁻¹Mpc (∼300 Mpc) on a side and contains 2500³ dark matter particles and an equal number of gas cells, implying a baryonic mass resolution of 7.44 × 10⁶ h⁻¹M_⊙ and a dark matter particle mass of 3.98 × 10⁷ h⁻¹M_⊙. Its outputs are publicly available at www.tng-project.org/data (Nelson et al. 2019). The differences in volume, resolution, and hydrodynamic implementation mean that the comparison between the two simulations can only be used as an approximate diagnostic and should not be interpreted as a full validation. Nevertheless, it provides a useful way to obtain a rough estimate of the amplitude of the bias introduced by the MTNG implementation.

In Fig. D.1, we compare the evolution of the gas mass as a function of m_sub/m_infall in the MTNG and TNG300 simulations, similarly to Figs. 3 and 5. We find that the gas masses associated with satellites in MTNG are systematically lower than those measured in TNG300, consistent with the expected bias introduced by the MTNG implementation of SUBFIND-HBT. However, the differences are not catastrophic and are of similar magnitude to the scatter shown in Fig. 3. In addition, the overall shape of the relation is broadly similar in the two simulations, and Eq. (1) is also capable of characterising the evolution of the gas mass in TNG300 (magenta dash-dotted line). We also examine the gas-mass profiles of satellite galaxies in TNG300 and find qualitatively similar behaviour to that seen in MTNG, in the sense that gas loss does not proceed strictly from the outside in (not shown here).

Fig. D.1. Gas-mass evolution in MTNG and TNG300. This figure is similar to the top right panel of Fig. 5, but shows only the galaxy sample with log(m200c/h−1M⊙) = 13.5 and includes the predictions from the TNG300 simulation. The symbols represent the data from the simulations, while the dash-dotted lines represent the analytic models described in Sect. 3.2.

These comparisons suggest that the main results of this work concerning the evolution of satellite gas content remain qualitatively robust, in the sense that the overall trend appears to be preserved. At the same time, it indicates that the absolute gas-mass scale inferred from MTNG should be interpreted with caution, as it is likely biased low for diffuse gas associated with satellite subhaloes.

A more detailed assessment of this effect, including a reanalysis of the main gas-related tests in TNG300 and/or with updated subhalo-identification procedures, is left for future work. The comparison presented here is intended only to provide an approximate estimate of the magnitude of the bias and to guide the interpretation of the gas-related results reported in this paper.

Appendix E: Best-fitting parameters and scatter of the analytic models

In this appendix, we present the best-fitting parameters for the analytic models describing the evolution of stellar mass as a function of the remaining maximum circular velocity fraction, gas mass as a function of the remaining subhalo mass fraction, and the r- and U-band magnitudes as a function of the remaining maximum circular velocity fraction. We also include a characterisation of the scatter of these relations and the root-mean-square residuals as a quantitative measure of the quality of the fits (Table E.1).

All Tables

All Figures

Fig. 1. Satellite HOD and stellar-mass relation. Top panel: The solid lines show the satellite HOD for galaxies in the MTNG simulation selected by peak subhalo mass, i.e. the maximum subhalo mass reached by each galaxy during its evolutionary history, while the dashed black line shows all subhaloes with peak subhalo mass above 1011 h−1 M⊙. The vertical coloured bands correspond to the halo-mass selection used throughout this work. Bottom panel: The solid line shows the median stellar mass of satellite galaxies as a function of peak subhalo mass, while the shaded regions represent the 16th and 84th percentiles. The vertical coloured bands correspond to the peak-subhalo-mass selection used throughout this work.

In the text

Fig. 2. Illustration of a satellite with mstell ∼ 1010.5 h−1 M⊙ infalling into a halo of ∼ 1013.5 h−1 M⊙ at z = 0. The solid circles show the satellite’s position relative to the halo centre. The size of each dot scales with the stellar mass, gas mass, star formation rate, and subhalo mass (in the first, second, fifth, and sixth columns, respectively), and with the difference in the r- and U-band magnitudes (in the third and fourth columns), relative to their values at the moment of infall. The larger dotted circle represents one virial radius (r200c).

In the text

Fig. 3. Evolution of the stellar mass, gas mass, and r- and U-band magnitudes of satellite galaxies at z = 0 respectively shown as a function of the remaining subhalo mass fraction (first column), pericentric distance (second column), time since infall (third column), and the remaining maximum circular velocity fraction (fourth column). The different colours represent different host halo masses, while the different line styles represent different peak subhalo masses, as labelled. For simplicity, we show different halo masses only for the intermediate-mpeak sample. The shaded region represents the 16th and 84th percentiles of the distribution.

In the text

	Fig. 4. Evolution of satellite galaxy properties at z = 1. This figure is similar to Fig. 3, but for z = 1.
In the text

Fig. 5. Evolution of the stellar mass as a function of vmax/vinfall (top left panel), the gas mass as a function of msub/minfall (top right panel), and the r- and U-band magnitudes as a function of vmax/vinfall (bottom left and right panels) for galaxies with 11.5 < log(mpeak/h−1 M⊙) < 12.0. The symbols represent the data from the MTNG simulation (as shown in Fig. 3), while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled.

In the text

Fig. 6. Mass profiles of satellite subhaloes with peak subhalo masses between 1011.5 and 1012 h−1 M⊙ (top row) and between 1012 and 1012.5 h−1 M⊙ (bottom row). The first, second, and third columns show the subhalo mass, gas mass, and stellar mass profiles, respectively. The solid coloured lines represent different levels of subhalo mass loss, while the dashed black lines represent the profiles of central subhaloes. Each profile is normalised by the comoving virial radius (r200c) of the satellite at the moment of infall.

In the text

	Fig. 7. Ratios of satellite mass profiles. This figure is similar to Fig. 6, but shows the ratio of the profiles to the 0.8 < msub/mpeak < 1 sample.
In the text

Fig. A.1. Gas fraction of central subhaloes in the MTNG simulation, normalised by Ωb/ΩM, as a function of the peak value of the maximum circular velocity reached by each galaxy during its evolutionary history (i.e. vpeak). The solid coloured line represents the median of the distribution, while the shaded region represents the 68th percentile. The dashed black line represents the fit to this distribution (see Appendix A for more details).

In the text

Fig. B.1. Stellar and gas masses of centrals and satellites shortly after infall. Top panel: The solid black line and the shaded regions correspond to the median of the distribution for central galaxies and the 16th and 84th percentiles, respectively. The orange line shows the median stellar mass for all satellite galaxies, while the blue line shows only satellites with msub/mpeak > 0.97 (representing satellites shortly after infall). Bottom panel: The same comparison is shown for gas mass normalised by peak subhalo mass. The lower panel in each figure shows the ratio of the coloured lines to the black line (central galaxies).

In the text

Fig. B.2. Halo bias in the MTNG-DMO simulation. The halo bias is computed as the ratio of the halo–matter cross-correlation function to the theoretical matter correlation function. The dotted and dashed lines show the biases of the haloes in the lower and upper 20% of the concentration distribution, respectively, in bins of halo mass (i.e. the concentration halo assembly bias). The shaded regions correspond to the bins in mpeak used in this work, which should approximately map onto similar halo-mass ranges.

In the text

Fig. B.3. Stellar and gas masses as a function of environment. Top panel: The solid black line and the shaded regions correspond to the median of the distribution for central galaxies and the 16th and 84th percentiles, respectively. The dashed and dotted lines represent central galaxies in haloes residing in the 20% most and least dense environments, respectively. The environmental classification is based on the individual bias-per-object and is performed in logarithmic bins of halo mass. Bottom panel: The same comparison is shown for gas mass normalised by peak subhalo mass.

In the text

	Fig. B.4. Mass profiles of centrals in dense and underdense environments. This figure is similar to Fig. 6, but for centrals in the 20% densest and least dense environments. The environmental classification is based on the individual bias-per-object and is performed in logarithmic bins of halo mass.
In the text

Fig. C.1. Evolution of stellar mass and magnitudes as a function of remaining subhalo mass fraction. This figure is similar to Fig. 5, but shows the evolution of the stellar mass and the r- and U-band magnitudes as a function of msub/minfall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0. The symbols represent the data from the MTNG simulation, as shown in Fig. 3, while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled.

In the text

Fig. C.2. Evolution of satellite properties as a function of time since infall. This figure is similar to Fig. 5, but shows the evolution of the stellar mass, gas mass, and r- and U-band magnitudes as a function of time since infall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0. The symbols represent the data from the MTNG simulation, as shown in Fig. 3, while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled.

In the text

Fig. C.3. Analytic fits to satellite galaxy evolution at z = 1. This figure is similar to Fig. 5, but shows the evolution of the stellar mass as a function of vmax/vinfall, the gas mass as a function of msub/minfall, and the r- and U-band magnitudes as a function of vmax/vinfall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0 at z = 1. The symbols represent the data from the MTNG simulation, as shown in Fig. 3, while the dash-dotted lines represent the analytic models described in Sect. 3.2. The different colours represent different host halo masses, as labelled. As a reference, the dotted orange line shows the evolution of the intermediate-halo mass sample at z = 0.

In the text

Fig. C.4. Evolution of the stellar mass as a function of vmax/vinfall, the gas mass as a function of msub/minfall, and the r- and U-band magnitudes as a function of vmax/vinfall for galaxies with 11.5 < log(mpeak/h−1M⊙) < 12.0. The different colours represent different host halo masses, as labelled. The symbols represent the data from the MTNG simulation, while the dash-dotted lines represent an analytic model fitted simultaneously to all halo masses.

In the text

	Fig. C.5. Evolution of stellar mass as a function of vmax/vinfall. The dash-dotted line represents an analytic model, obtained through symbolic regression, fitted simultaneously to all halo masses.
In the text

	Fig. D.1. Gas-mass evolution in MTNG and TNG300. This figure is similar to the top right panel of Fig. 5, but shows only the galaxy sample with log(m200c/h−1M⊙) = 13.5 and includes the predictions from the TNG300 simulation. The symbols represent the data from the simulations, while the dash-dotted lines represent the analytic models described in Sect. 3.2.
In the text