Orbital dynamics and histories of satellite galaxies around Milky Way – mass galaxies in the FIRE simulations

ABSTRACT

1 INTRODUCTION

The satellite galaxies around the Milky Way (MW) and M31 in the Local Group (LG) are unique systems in that we can measure both resolved stellar populations and full orbits to derive orbital histories. Key answerable questions regarding satellite formation histories include: What are their orbits today and what were their orbital histories? When did they first become satellites? How many were satellites of a lower-mass group, like the Large Magellanic Cloud (LMC; Kallivayalil et al. 2018; Pardy et al. 2020; Patel et al. 2020), and when did they become satellites of such groups? How close have they orbited to the MW and M31? Understanding the answers to these questions within the LG will also improve our understanding of satellite evolution in systems beyond the LG.

Thanks to numerous studies (e.g. Kallivayalil et al. 2013, 2018; Fritz et al. 2018b; Gaia Collaboration 2018; Patel et al. 2020) and HST treasury programs (e.g. GO-14734, PI Kallivayalil; GO-15902, PI Weisz), we now have or soon will have 3D velocity information for the majority of the satellites in the LG, with continued higher precision (e.g. see updated values between McConnachie 2012; Fritz et al. 2018a, and references therein). With the advent of new data coming from studies such as the Satellites Around Galactic Analogs (SAGA) survey (Geha et al. 2017; Mao et al. 2021), which focus on measuring properties of satellites of MW-mass galaxies beyond the LG, understanding general orbital/infall histories of satellites is imperative. Given this rich data, combined with observationally informed estimates of the MW’s gravitational potential (e.g. Bland-Hawthorn & Gerhard 2016; McMillan 2017; Li et al. 2020a; Deason et al. 2021; Correa Magnus & Vasiliev 2022), one can constrain the orbital histories of satellites. However, a key challenge is understanding limitations and biases from assuming a static, non-evolving potential (e.g. Klypin et al. 1999; D’Souza & Bell 2022).

One of the most important events in a satellite galaxy’s history is when it first orbited within the MW-mass halo, or any more massive halo. After falling into a more massive host halo, a satellite galaxy temporarily can orbit outside of the host’s virial radius (a ‘splashback’ satellite); typically such satellites remain gravitationally bound and orbit back within the host’s virial radius (e.g. Ludlow et al. 2009; Wetzel et al. 2014). As a satellite orbits within the host halo, its hot gas can quench the satellite’s star formation through ram-pressure stripping of gas (e.g. Gunn & Gott 1972; van den Bosch et al. 2008; Fillingham et al. 2019; Rodriguez Wimberly et al. 2019; Samuel et al. 2022). Not only can a lower-mass galaxy fall into a MW-mass halo, but it can become a satellite of another intermediate-mass galaxy before falling into the MW-mass halo, called ‘pre-processing’ (e.g. D’Onghia & Lake 2008; Li & Helmi 2008; Deason et al. 2015; Wetzel, Deason & Garrison-Kimmel 2015; Kallivayalil et al. 2018; Jahn et al. 2019; Patel et al. 2020). This process can suppress and even quench star formation in the low-mass galaxy before it falls into its MW-mass halo (e.g. Jahn et al. 2022; Samuel et al. 2022). If satellites fell into a MW-mass halo via a group, they should have similar orbits, at least for one or a few orbital time-scales, with broadly similar orbital angular momenta and energy (e.g. Sales et al. 2011; Deason et al. 2015; Sales et al. 2017; Pardy et al. 2020). Some theoretical studies suggest that no current satellites of the MW were satellites during the epoch of reionization (e.g. Wetzel et al. 2015; Rodriguez Wimberly et al. 2019) at z ≳ 6; thus, if the satellites quenched at z ≳ 6, the host environment could not have quenched them, such that the effects of the host environment and cosmic reionization are separable, in principle.

Many works have studied infall histories and orbital properties of simulated satellite galaxies of MW-mass haloes (e.g. Slater & Bell 2013; Wetzel et al. 2015; Li et al. 2020b; Bakels, Ludlow & Power 2021; D’Souza & Bell 2021; Robles & Bullock 2021; Ogiya, Taylor & Hudson 2021), sometimes with the intent of deriving properties of observed satellites of the MW (e.g. Rocha, Peter & Bullock 2012; Fillingham et al. 2019; Miyoshi & Chiba 2020; Rodriguez Wimberly et al. 2022). However, many such previous works used dark-matter-only (DMO) simulations, and the inclusion of baryonic physics is critical to model accurately the satellite population (see Brooks & Zolotov 2014; Bullock & Boylan-Kolchin 2017; Sales, Wetzel & Fattahi 2022).

Because dynamical friction robs (primarily higher-mass) satellites of their orbital energy, one might expect satellites to shrink monotonically in orbital radius over time, until the satellite merges/disrupts. Many studies implement idealized simulations (both with and without baryonic physics) that incorporate the effects of dynamical friction, mass-loss, and ram-pressure stripping (e.g. Weinberg 1986; Taylor & Babul 2001; Peñarrubia, Kroupa & Boily 2002; Peñarrubia & Benson 2005; Amorisco 2017; Jiang et al. 2021), but testing this assumption in a cosmological setting is imperative.

Additionally, because the LMC has satellites of its own (e.g. D’Onghia & Lake 2008; Deason et al. 2015; Kallivayalil et al. 2018), many studies have investigated satellite spatial distributions and their dynamics with an additional LMC-like contribution to the host potential to test the dynamical effects of a massive satellite on nearby lower-mass galaxies (e.g. Garavito-Camargo et al. 2019; Patel et al. 2020; Samuel et al. 2021; Vasiliev, Belokurov & Erkal 2021; D’Souza & Bell 2022; Pace, Erkal & Li 2022).

In this paper, we examine the orbital histories of satellite galaxies using baryonic simulations that match general observed properties of satellites of MW-mass galaxies. Specifically, we use the FIRE-2 cosmological zoom-in simulations, which form realistic MW-mass galaxies (e.g. Garrison-Kimmel et al. 2018; Hopkins et al. 2018; Sanderson et al. 2020; Bellardini et al. 2021) and populations of satellite galaxies around them (Wetzel et al. 2016; Garrison-Kimmel et al. 2019a; Samuel et al. 2020, 2022). The two main topics that we explore are: (1) The relation of orbital properties of satellite galaxies at z = 0 to their orbital histories, including lookback times of infall, distances from the MW-mass host, and stellar masses. These relations not only help characterize the orbits of satellites, but their orbital histories also provide insight, or caution, into approximating history-based properties, such as infall time or pericentre information, based on their present-day properties, such as distance and stellar mass. (2) Testing a common expectation that the orbits of satellite galaxies shrink over time, that is, that a satellite’s most recent pericentric distance is the minimum that it has experienced.

In Santistevan et al. (in preparation), we will compare directly the orbits of satellites in cosmological simulations of MW-mass galaxies to orbits from integration in a static, idealized MW-mass halo (see also Vasiliev et al. 2021; D’Souza & Bell 2022, and references therein).

2 METHODS

2.1 FIRE-2 simulations

We use the cosmological zoom-in baryonic simulations of both isolated MW/M31-mass galaxies and LG-like pairs from the Feedback In Realistic Environments (FIRE) project¹ (Hopkins et al. 2018). We ran these simulations using the hydrodynamic plus N-body code gizmo (Hopkins 2015), with the mesh-free finite-mass hydrodynamics method (Hopkins 2015). We used the FIRE-2 physics model (Hopkins et al. 2018) that includes several radiative heating and cooling processes such as Compton scattering, bremsstrahlung emission, photoionization and recombination, photoelectric, metal-line, molecular, fine-structure, dust-collisional, and cosmic ray heating across temperatures 10–10¹⁰ K, including the spatially uniform and redshift-dependent cosmic ultraviolet (UV) background from Faucher-Giguère et al. (2009), for which H i reionization occurs at z_reion ≈ 10. Star formation occurs in gas that is self-gravitating, Jeans unstable, molecular (following Krumholz & Gnedin 2011), and dense (n_H > 1000 cm⁻³). Star particles represent single stellar populations, assuming a Kroupa (2001) initial mass function, and they evolve along stellar population models from starburst99v7.0 (Leitherer et al. 1999), inheriting masses and elemental abundances from their progenitor gas cells. FIRE-2 simulations also include the following stellar feedback processes: core-collapse and Type Ia supernovae, stellar winds, and radiation pressure.

We used the code music (Hahn & Abel 2011) to generate the cosmological zoom-in initial conditions at z ≈ 99 within periodic cosmological boxes of length 70.4–172 Mpc, sufficiently large to avoid unrealistic periodic gravity effects on individual MW-mass haloes. For each simulation, we save 600 snapshots with 20–25 Myr spacing down to z = 0, assuming a flat ΛCDM cosmology. Consistent with Planck Collaboration (2020), we used cosmological parameters in the following ranges: h = 0.68–0.71, σ₈ = 0.801–0.82, n_s = 0.961–0.97, Ω_Λ = 0.69–0.734, Ω_m = 0.266–0.31, and Ω_b = 0.0449–0.048.

Our galaxy sample consists of the 12 MW/M31-mass galaxies in Santistevan et al. (2020), as well as one additional galaxy, ‘m12r’, first introduced in Samuel et al. (2020). These are from the Latte suite of isolated MW/M31-mass galaxies introduced in Wetzel et al. (2016) and the ‘ELVIS on FIRE’ suite of LG-like MW+M31 pairs, introduced in Garrison-Kimmel et al. (2019a). Table 1 lists their stellar mass, M_star,90, halo mass, M_200m, and radius R_200m at z = 0, where both are defined at 200 times the matter density at z = 0. The Latte suite consists of haloes with |$M_{\rm 200m} = 1\mathrm{ - }2 \times 10^{12} \, {\rm M}_{\odot }$| at z = 0 with no other similar-mass haloes within 5 × R_200m. We also chose m12r and m12w to have LMC-mass satellite analogues near z ≈ 0 (Samuel et al. 2020). The initial masses of star particles and gas cells is |$7100 \, {\rm M}_{\odot }$|⁠, but the average star particle mass at z = 0 is |$\approx 5000 \, {\rm M}_{\odot }$| from stellar mass-loss. Within the zoom-in region, the mass of DM particles is |$3.5 \times 10^4 \, {\rm M}_{\odot }$|⁠. The gravitational softening lengths are 4 and 40 pc (Plummer equivalent), co-moving at z > 9 and physical thereafter, for star and DM particles, respectively. Gas cells use adaptive force softening, equal to their hydrodynamic smoothing, down to 1 pc.

Each pair of haloes in the ‘ELVIS on FIRE’ suite of LG-like pairs was chosen based on their individual mass (⁠|$M_{\rm 200m} = 1 \mathrm{ -} 3 \times 10^{12} \, {\rm M}_{\odot }$|⁠) and combined masses (total LG mass between |$2 \mathrm{ -} 5 \times 10^{12} \, {\rm M}_{\odot }$|⁠), as well as their current separation (⁠|$600 \mathrm{ - }1000\, {\rm kpc}$|⁠) and radial velocities at z = 0 (⁠|$\rm \upsilon _{rad} \lt 0$|⁠), and isolated environment (no other massive haloes within |$2.8\, {\rm Mpc}$| of either host centre). The mass resolution is ≈2× better in the ‘ELVIS on FIRE’ suite, with initial baryonic particle masses |$3500 \mathrm{ -} 4000 \, {\rm M}_{\odot }$|⁠. For all results in this paper, we investigated possible differences between the isolated and LG-like MW-mass galaxies, which also partially tests resolution convergence, and we find negligible differences between the two samples.

These 13 host galaxies reflect formation histories of general MW/M31-mass (or LG-like) galaxies within our selection criteria and exhibit observational properties broadly similar to the MW and M31, including: realistic stellar haloes (Bonaca et al. 2017; Sanderson et al. 2018), dynamics of metal-poor stars from early galaxy mergers (Santistevan et al. 2021), satellite galaxy stellar masses and internal velocity dispersions (Wetzel et al. 2016; Garrison-Kimmel et al. 2019b), radial and 3D spatial distributions (Samuel et al. 2020, 2021), and star-formation histories and quiescent fractions (Garrison-Kimmel et al. 2019b; Samuel et al. 2022).

2.2 Halo/Galaxy catalogues and merger trees

We use the rockstar 6D halo finder (Behroozi, Wechsler & Wu 2013a) to generate (sub)halo catalogues using only DM particles at each of the 600 snapshots, and CONSISTENT-TREES (Behroozi et al. 2013b) to generate merger trees. As a consequence of the large zoom-in volume for each host, there is no low-resolution DM particle contamination in any of the (sub)haloes that we analyse.

Samuel et al. (2020) describe our star particle assignment to (sub)haloes in post-processing; we briefly review it here. We first select star particles within d < 0.8R_halo (out to a maximum distance of |$30 \, {\rm kpc}$|⁠) with velocities υ < 2V_circ,max of the (sub)halo’s centre-of-mass (COM) velocity. We then keep only the star particles within d < 1.5R_star,90 (the radius enclosing 90 per cent of the stellar mass) of the (then) current member stellar population’s COM and halo centre position. We further kinematically select the star particles with velocities υ < 2σ_vel,star (the velocity dispersion of current member star particles) of the COM velocity of member star particles. Finally, we iterate on these spatial and kinematic criteria, which guarantees that the COM of the galaxy and (sub)halo are consistent with one another, until the (sub)halo’s stellar mass converges to within 1 per cent.

We use two publicly available analysis packages: HaloAnalysis² (Wetzel & Garrison-Kimmel 2020a) for assigning star particles to haloes and for reading and analysing halo catalogues/trees, and GizmoAnalysis³ (Wetzel & Garrison-Kimmel 2020b) for reading and analysing particles from gizmo snapshots.

2.3 Selection of satellites

We include all luminous satellite galaxies at z = 0 with |$M_{\rm star} \gt 3\times 10^4\, {\rm M}_{\odot }$| that have crossed within their MW-mass host halo’s R_200m(z). This lower limit on stellar mass corresponds to ≈6 star particles, the limit for reasonably resolving the total stellar mass (Hopkins et al. 2018). Our sample includes ‘splashback’ satellites that are currently beyond the host’s R_200m, which are typically still gravitationally bound to the host but simply near apocentre (e.g. Wetzel et al. 2014). As Table 1 shows, the number of surviving luminous satellites at z = 0, including this splashback population, per host ranges from 26 to 52, and our sample totals 473 satellites. Both Garrison-Kimmel et al. (2014) and Wetzel et al. (2015) showed that in the ELVIS DMO simulation suite, the average number of subhaloes that would typically host galaxies with |$M_{\rm star}\gtrsim 10^5\, {\rm M}_{\odot }$| is ∼31–45. However, because stellar feedback and baryonic physics can affect galaxy formation, Wetzel et al. (2016) and Garrison-Kimmel et al. (2019a) showed that the number of satellites above this mass range decreases to ∼13–15. More recently, Samuel et al. (2020) showed that the radial distributions of these satellites are consistent with the MW and M31 out to |$300\, {\rm kpc}$|⁠. The MW and M31 each have 13 and 27 satellites, respectively, and the MW-mass hosts in our simulations bracket these values with 11 to 27 satellites. Unless otherwise stated, in our analysis we refer to luminous satellites, i.e. satellites containing stars, as simply ‘satellites’.

In computing host-averaged results below, to avoid biasing our results to the hosts with larger satellite populations, we oversample the satellites so that each host contributes a nearly equal fraction of satellites to the total. Specifically, we multiply the number of satellites in the MW-mass host with the largest population (Romulus, with 52 satellites) by 10, which results in an oversampled population of 520. Then, for each of the other MW-mass hosts, we divide 520 by the number of their satellites and obtain the nearest integer multiplicative factor, m, that we apply to each host’s satellite population, N_sat (see Table 1), so that each host contains ≈500–530 satellites, or that their satellite populations are within 5 per cent of one another. Thus, when plotting properties, such as pericentre distances, we count each satellite in a given host m times for each property in the figures.

Fig. 1 shows the relation between stellar mass and halo mass (SMHM) for our satellites. We show stellar mass atz= 0 versus peak (sub)halo mass throughout its history. We include the median SMHM relation of non-satellites in the same simulations, which we define as low-mass galaxies that never crossed the virial radius of the MW-mass host, and that currently orbit beyond |$1\, \rm {Mpc}$| at z = 0, in the blue dot–dashed line. The formation histories of non-satellites differ from satellites because they form in less dense regions of the Universe. This, along with the UV heating of gas in isolated low-mass galaxies, may explain their slightly smaller M_star; however, a deeper investigation is outside of the scope of this paper. The grey dotted line at |$M_{\rm star} = 3\times 10^4\, {\rm M}_{\odot }$| shows our lower limit in stellar mass.

Figure 1.

Stellar mass, M_star, at z = 0 versus peak halo mass, |${M}_\rm {halo,peak}$|⁠, for all 473 satellites with |$M_{\rm star} \gt 3 \times 10^4 \, {\rm M}_{\odot }$| across all 13 host galaxies. The solid blue line shows the median, and the dark and light shaded regions show the 68th percentile and full distribution, respectively. We compare this with the stellar-halo mass relation of non-satellite low-mass galaxies that never fell into the MW-mass host and currently orbit beyond |$1\, \rm {Mpc}$| at z = 0 (blue dot–dashed). The dotted grey line indicates the minimum M_star in our sample. For comparison, we also show extrapolations of the stellar-halo mass relations from Moster, Naab & White (2013) (red), Garrison-Kimmel et al. (2017a) (green), and Behroozi et al. (2020) (black), which broadly agree with our sample at |${M}_\rm {halo,peak}\gtrsim 10^9\, {\rm M}_{\odot }$|⁠. Additionally, we show the 68th percentile and full distribution of values from higher resolution, isolated low-mass galaxies from Fitts et al. (2017) in the dark and light-shaded pink regions, respectively, and the median across their sample in the horizontal pink line. We finally show four higher resolution, isolated low-mass galaxies from Wheeler et al. (2019) as stars; however, the galaxy with the smallest M_star in this sample is beyond our resolution limit. At all halo masses, the 68th percentile in satellite M_star is 0.5–1 dex, with a smaller range at the extreme masses from lower statistics. The relation flattens at |$M_{\rm star} \lesssim 10^5 \, {\rm M}_{\odot }$| simply because of our stellar mass limit; therefore, our fiducial sample is complete in M_star to |$\gt 3 \times 10^4 \, {\rm M}_{\odot }$|⁠, while we are complete in |${M}_\rm {halo,peak}$| to |$\gtrsim 3 \times 10^9 \, {\rm M}_{\odot }$| (see Appendix  A for results for selecting satellites via |${M}_\rm {halo,peak}$|⁠).

Open in new tabDownload slide

To compare to other SMHM relations from the literature, we also show extrapolations from Moster et al. (2013), Garrison-Kimmel et al. (2017a), and Behroozi et al. (2020) in the dashed red, green, and black lines, respectively. We also include the median M_star, 68th percentile, and full distribution of values from Fitts et al. (2017) of higher resolution (initial baryon masses of |$500\, {\rm M}_{\odot }$|⁠), isolated low-mass galaxies in the pink horizontal line, and dark- and light-pink-shaded regions, respectively. Finally, we include four higher resolution, isolated low-mass galaxies from Wheeler et al. (2019; named m09, m10q, m10v, and m10vB in their work), with initial baryon masses of |$30\, {\rm M}_{\odot }$| (magenta stars). The SMHM relation in our sample broadly agrees with these (extrapolated) semi-empirical estimates and values of isolated low-mass galaxies; for a more detailed discussion about the SMHM relation in FIRE-2, see Hopkins et al. (2018) and Wheeler et al. (2019). The low-mass end of our SMHM relation in Fig. 1 flattens at |$M_{\rm star}\lesssim 10^5\, {\rm M}_{\odot }$|⁠, purely because of our stellar mass selection of |$M_{\rm star}\gt 3\times 10^4\, {\rm M}_{\odot }$|⁠. We note that the galaxy with the smallest M_star from Wheeler et al. (2019) is beyond the resolution limit in our sample, and the second smallest M_star galaxy would be only marginally resolved. However, with better resolution, the SMHM relation in Fig. 1 would likely follow a similar trend as the extrapolations and it is likely that the lowest M_star isolated galaxy would lie in the lower end of the full distribution scatter. Thus, while our sample is complete in stellar mass, we are complete in halo mass for |${M}_\rm {halo,peak}\gtrsim 3\times 10^9\, {\rm M}_{\odot }$|⁠. We find only minor differences in our results if selecting satellites via |${M}_\rm {halo,peak}$| instead (see Appendix  A).

We checked the SMHM relation in Fig. 1 instead using the peak stellar mass throughout a galaxy’s history, M_star,peak. The two relations are similar, but the SMHM relation with M_star,peak has a |$\mathrm{ \approx} 5{{\ \rm per\ cent}}$| higher normalization, on average, because of stellar mass-loss after infall.

2.4 Numerical disruption

Many previous studies (e.g. van den Bosch & Ogiya 2018; van den Bosch et al. 2018; Webb & Bovy 2020) noted that without proper mass and spatial resolution, satellite subhaloes can suffer from artificial numerical disruption too quickly. Thus, sufficient resolution and implementation of the relevant physics is necessary to model accurately the evolution of satellite galaxies. As a partial test of this, we investigated differences between the isolated and LG-like satellite populations, which have a ≈2 × resolution difference, and we saw no strong differences in our results. Samuel et al. (2020) also tested resolution convergence of the satellite radial profiles using the FIRE-2 simulations with dark matter particle masses of |$m_{\rm dm} = 3.5 \times 10^4\, {\rm M}_{\odot }$| and |$m_{\rm dm} = 2.8 \times 10^5\, {\rm M}_{\odot }$| and found generally consistent results between the two, because there are enough DM particles (≳ 2 × 10⁴) in the lowest-mass luminous subhaloes (⁠|${M}_\rm {halo,peak}\gtrsim 10^8 \, {\rm M}_{\odot }$|⁠) to prevent numerical disruption, and many more than that in our typical (more massive) luminous subhaloes, thus satisfying criteria such as in van den Bosch & Ogiya (2018). Also, our DM particle mass resolution is |$m_{\rm dm} = 2 \mathrm{ -} 3.5\times 10^4\, {\rm M}_{\odot }$| and DM force softening length is 40 pc, significantly better than previous work like Wetzel et al. (2015), who used the ELVIS DMO simulations with |$m_{\rm dm} = 1.9 \times 10^5 \, {\rm M}_{\odot }$| and 140 pc but found results broadly consistent with ours. Nonetheless, any simulations like ours necessarily operate at finite resolution, which inevitably leads to some degree of numerical disruption, which any reader should bear in mind.

2.5 Calculating pericentre

We calculate pericentres by tracking the main progenitor back in time using the merger trees, which store each satellite’s galactocentric distance at each of the 600 snapshots. We first ensure that the satellite is within the MW-mass host halo at a given snapshot. Then, we check if the galactocentric distance reaches a minimum within ±20 snapshots, corresponding to a time window of |$\mathrm{ \approx }1 \, {\rm Gyr}$|⁠. Given the |$\mathrm{ \approx} 25\, {\rm Myr}$| time spacing between snapshots, we fit a cubic spline to the distance and time arrays across this interval. We then find the minimum in the spline-interpolated distance(s), and record the corresponding spline-interpolated time.

We checked how our results differ if varying the window to ±4, 8, and 10 snapshots by visually inspecting each satellite’s orbit history and conclude that a window of ±20 snapshots reduces nearly all ‘false’ pericentres, that is, instances in which the criteria above are met because of numerical noise in the orbit or a short-lived perturbation. Because the COM of the MW-mass galaxy does not perfectly coincide with the COM of its DM host halo, we also checked how our results vary between using distances with respect to the centre of the galaxy versus the halo: we find no meaningful difference. We additionally checked how our results vary when using distances with respect to the centre of the satellite galaxy versus the centre of the satellite (sub)halo, again finding no significant differences.

2.6 Calculating gravitational potential and total energy

We also explore trends with a satellite’s specific orbital total energy, E at z = 0, defined as the sum of the kinetic and potential energies. Our simulations store the value of the gravitational potential for all particles at each snapshot, so to calculate the potential for each satellite at z = 0, we select all star, gas, and DM particles within |$\pm 5 \, {\rm kpc}$| of the satellite (sub)halo’s virial radius, to limit biasing from the satellite’s self-potential, and we compute the mean potential across these particles. Given that some satellites are in LG-like environments, we normalize E at the MW-mass halo radius, such that E(d = R_200m) = 0, that is, the sum of the host potential at R_200m and the kinetic energy of a circular orbit at R_200m is 0.

3 RESULTS

Throughout the paper, we present results for all satellite galaxies at z = 0, based on their stellar mass (⁠|$M_{\rm star} \gt 3 \times 10^4 \, {\rm M}_{\odot }$|⁠), across all of our MW-mass hosts. Although Santistevan et al. (2020) noted that MW-mass haloes/galaxies in LG-like pairs formed |$\mathrm{ \sim }1.6 \, {\rm Gyr}$| earlier than those that are isolated, we compared satellites based on isolated versus LG-like environment and find negligible differences in any properties that we investigate. This agrees with the lack of dependence in Wetzel et al. (2015), who investigated the infall histories of satellite subhaloes in the ELVIS DMO simulations. Appendix  A examines how our results change by selecting satellites by their peak halo mass. In summary, we find qualitatively similar results for selecting via stellar mass, but given the scatter in the SMHM relation, the trends with halo mass are smoother and any features are sharper. Although the DM (sub)halo mass is more dynamically relevant to the orbits of satellite galaxies, we present our results with a stellar-mass selected sample, because it is observationally easier to measure than halo mass. Finally, Appendix  B compares our results for baryonic simulations against DMO simulations of the same systems. In summary, the lack of a MW-mass galaxy in the DMO simulations allows satellites to survive longer, orbit closer to the centre of the halo, and complete more orbits.

We examine trends guided by which orbital properties are relevant to different phenomena. As we will show, specific angular momentum and specific total energy provide insight into when a satellite fell into the MW-mass halo. We also explore trends with satellite mass, in part to understand where dynamical friction becomes important: from equation (1), for the MW-mass haloes with |$M_{\rm 200m}(z = 0) \approx 10^{12}\, {\rm M}_{\odot }$| (and lower M_200m at earlier times), we expect dynamical friction to significantly affect satellites with |$M_{\rm 200m} \gtrsim 3 \times 10^{10}\, {\rm M}_{\odot }$|⁠, or |$M_{\rm star} \gtrsim 10^8 \, {\rm M}_{\odot }$|⁠. We also focus on infall times, to understand how long satellites have been orbiting in the host halo environment, and we explore the incidence of pre-processing in a lower-mass group prior to MW-mass infall. We also examine properties of orbital pericentre, given that satellites typically feel the strongest gravitational tidal force and the strongest ram pressure at pericentre.

In this paper, we present trends for the simulated satellite populations only; however, in the future, we plan to investigate differences in the simulations and results obtained from idealized orbit modelling methods. Ultimately, we will provide a framework to derive similar orbital properties from satellites in the MW and M31 using the satellite populations in the simulations, and compare the results. Thus, we leave direct observational comparisons for future work.

3.1 Orbital properties today

We first investigate the instantaneous orbital properties of satellites at z = 0, including approximate integrals of motion like angular momentum and energy. Fig. 2 shows total velocity, specific orbital total energy, E, and specific angular momentum, ℓ, as a function of lookback time of infall into the MW-mass halo, |$t^{\rm lb}_{\rm infall,MW}$|⁠, galactocentric distance, d, and stellar mass, M_star. The top of each column shows distributions of |$t^{\rm lb}_{\rm infall,MW}$|⁠, d, and M_star. In particular, the distribution of |$t^{\rm lb}_{\rm infall,MW}$| is relatively flat, with a modest peak |$\mathrm{ \approx} 9 \, {\rm Gyr}$| ago. For reference, all panels versus |$t^{\rm lb}_{\rm infall,MW}$| (left) include a vertical shaded region to represent the free-fall time at R_200m, |$t_{\rm ff} = \sqrt{ 3 \pi / \left(32 G \rho _{\rm 200m}\right) }$|⁠, where |$t_{\rm ff} = 2.8 \mathrm{ -} 3 \, {\rm Gyr}$| across our hosts. In all panels versus d (middle), the vertical shaded region shows the range of R_200m for the MW-mass haloes, as Table 1 lists. Similarly, all panels versus total velocity and ℓ show horizontal shaded bands that represent the range of V_200m and L_200m across hosts, where |$V_{\rm 200m} = \sqrt{GM_{\rm 200m}/R_{\rm 200m}}$| is the velocity of a circular orbit at the virial radius, and L_200m = V_200m × R_200m is the specific angular momentum of that orbit. Because the satellites fell in at different times and orbit at various distances with unique stellar masses, we do not expect them to have values equal to V_200m or L_200m, so we provide these shaded regions as a reference only.

Figure 2.

Instantaneous orbital dynamics of satellite galaxies at z = 0 versus their lookback time of infall into the MW-mass halo (left), distance from MW-mass host, d (middle), and stellar mass, M_star (right). The solid lines show the median for all satellites across our 13 MW-mass hosts, and the dark and light shaded regions show the 68th percentile and full distribution, respectively. At the top of each column, we show the histogram of these properties for the sample. We show common reference points to these properties for typical host haloes within our mass range in the horizontal and vertical grey shaded bands in some of these panels, such as the free-fall time, and virial radii, velocities, and angular momenta; see text for details. Top row: orbital total velocity. The median velocity for satellites that fell in |$\gt 4 \, {\rm Gyr}$| ago increases with earlier infall time from |$100\mathrm{-}200 \, \rm {km}\, \rm {s}^{-1}$| (top left). However, more recently infalling satellites show a peak near |$\approx 1.5 \, {\rm Gyr}$| of |$\approx 230 \, \rm {km}\, \rm {s}^{-1}$|⁠, highlighting satellites that are nearing their first pericentre, before decreasing again to |$\approx 100 \, \rm {km}\, \rm {s}^{-1}$| for satellites that only recently fell in. Velocity also decreases with host distance and is relatively flat with M_star at roughly |$135\, \rm {km}\, \rm {s}^{-1}$|⁠. Middle row: specific orbital total energy, E. We normalize E at the host’s R_200m so E(d = R_200m) = 0. E decreases for earlier-infalling satellites, with median values ranging from 0 to |$\mathrm{ -}3.8\times 10^4\, \rm {km}^{2}\, \rm {s}^{-2}$|⁠. Satellites that fell in earlier orbit deeper in the host halo gravitational potential and therefore are more bound (middle left). Satellites with |$M_{\rm star} \gtrsim 10^8 \, {\rm M}_{\odot }$| are more bound, because of the stronger dynamical friction they experience, and satellites closer to the host are more bound, because they orbit deeper in the gravitational potential (middle right). Bottom row: orbital specific angular momentum, ℓ. Specific angular momentum decreases for earlier-infalling satellites, from ≈3 to |$0.9 \times 10^4 \, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$| (bottom left), because satellites that fell in earlier fell in (and orbit) at smaller distances (see Fig. 5 and bottom left panel of Fig. 6). ℓ necessarily increases with d (bottom middle) and depends only weakly on satellite M_star (bottom right).

Open in new tabDownload slide

3.1.1 Total velocity

We first present trends in the total velocity, corresponding to the top row in Fig. 2. Considering the trends in total velocity with infall time (top left), satellites that fell in |$\lt 1 \, {\rm Gyr}$| ago have not yet experienced a pericentre, so they show an increase in total velocity with time since infall from |$0.5 \mathrm{ -} 1.5 \, {\rm Gyr}$| ago, because these satellites are near pericentre. The total velocity then decreases with increasing time since infall, because those satellites are now near apocentre. We see a similar, but weaker, peak in total velocity |$\mathrm{ \sim }6 \, {\rm Gyr}$| ago, from a marginally phase-coherent population near pericentre, but after a few orbits, satellites become out of phase. Satellites that fell in |$\lesssim 3 \, {\rm Gyr}$| ago typically have total velocities of |$150 \mathrm{ -} 250 \, \rm {km}\, \rm {s}^{-1}$|⁠, while earlier infalling satellites are only orbiting at |$100 \mathrm{ -} 185 \, \rm {km}\, \rm {s}^{-1}$|⁠. Comparing these infall times to t_ff in the vertical shaded bands, satellites with |$t^{\rm lb}_{\rm infall,MW}\lt t_{\rm ff}$| often have larger velocities, and again are likely to be on first infall, because they have not had enough time to reach pericentre. For reference, we also compare the satellite total velocities to the host’s virial velocity, V_200m (horizontal grey shaded band). Satellites that fell in |$t^{\rm lb}_{\rm infall,MW}\approx 3 \mathrm{ -} 11 \, {\rm Gyr}$| ago typically have comparable total velocities to the virial velocity, but the full scatter extends to both much larger and smaller values.

Next, the median velocity decreases with increasing d (top middle), from as high as |$260\, \rm {km}\, \rm {s}^{-1}$| for the closest satellites to |$80\, \rm {km}\, \rm {s}^{-1}$| for satellites near R_200m. The shaded region of the 68th percentile follows the median trend and the width is roughly constant at |$\approx 95 \, \rm {km}\, \rm {s}^{-1}$| across all distances. At R_200m, the total velocities of satellites are typically lower than V_200m both because they are not on perfectly circular orbits and because the population at large d is likely biased to lower values from the splashback population that typically have negligible velocities.

The median is nearly constant with M_star (top right), ranging from |$120 \,\,\mathrm{ to}\,\, 160 \, \rm {km}\, \rm {s}^{-1}$| for |$M_{\rm star} \approx 10^{4.75-8.25} \, {\rm M}_{\odot }$|⁠, which decreases to |$70 \mathrm{ - }100 \, \rm {km}\, \rm {s}^{-1}$| at higher mass, likely because sufficiently massive satellites experience significant dynamical friction that slows their orbit. Across all stellar masses, the average total velocity is |$\approx 135 \, \rm {km}\, \rm {s}^{-1}$| and the typical range of the 68th percentile is |$90 \mathrm{ - }205 \, \rm {km}\, \rm {s}^{-1}$|⁠. The median for all satellites with |$M_{\rm star} \lesssim 10^{8.5}\, {\rm M}_{\odot }$| show consistent values with V_200m.

3.1.2 Specific total energy

Next, the middle row in Fig. 2 shows trends in the specific orbital energy, E. Given that some satellites are in LG-like environments, which complicates computing the specific total energy beyond R_200m, we normalize E at the MW-mass halo’s R_200m so E(d = R_200m) = 0. Thus, satellites with E > 0 are the splashback population with apocentres beyond R_200m and are essentially all bound to the host halo.

The middle left panel shows that earlier infalling satellites are on more bound orbits, with the median E decreasing from ≈0 to |$-3.9 \times 10^4\, \rm {km}^{2}\, \rm {s}^{-2}$|⁠. This reflects the growth of the MW-mass halo over time. E increases with d from the host, so satellites are more bound at smaller distances. Thus, the median E is not constant with d, but this does not imply that specific energy is not being conserved over an orbit. As we explore below, |$t^{\rm lb}_{\rm infall,MW}$| correlates with d, such that satellites at smaller d fell in earlier (though with large scatter; see Fig. 5). This is largely because they fell in when the host R_200m was smaller, so they necessarily orbit at smaller d. This then leads to the correlation of E with d across the population at z = 0 (Fig. 2, centre panel).

Similar to the trends in total velocity, E does not strongly depend on M_star, except at |$M_{\rm star}\gtrsim 10^8\, {\rm M}_{\odot }$|⁠, where satellites experience significant dynamical friction, causing their orbits to become more bound, despite the fact that higher-mass satellites fell in more recently, as we show below.

3.1.3 Specific angular momentum

Last, we present trends in specific angular momentum in the bottom row of Fig. 2. The median specific angular momentum, ℓ (bottom), decreases across time since infall from ≈3 to |$1 \times 10^4 \, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$| (bottom left panel). Between |$t^{\rm lb}_{\rm infall,MW}\approx 2 \mathrm{ -} 6 \, {\rm Gyr}$|⁠, ℓ is nearly constant at |$\approx 2 \times 10^4 \, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$|⁠, which as we will show in Fig. 6 corresponds to satellites that completed 1–2 pericentres. The median ℓ is much higher for satellites that are to the left of the t_ff band, consistent with the higher velocities in these satellites. Similar to the top row, we define the reference host halo angular momentum, L_200m = V_200m × R_200m, and we show the range of values across hosts in the grey horizontal band. As expected, essentially all satellites have ℓ < L_200m, and most have ℓ significantly lower.

Generally, ℓ and its scatter increase with increasing d, as expected, given that ℓ = υ_tand, where υ_tan is the satellite’s tangential velocity. Because the MW-mass halo grew over time, satellites fell into their MW-mass halo on larger orbits at later times, which explains the increasing ℓ with more recent infall times.

The median ℓ is nearly constant with stellar mass at |$M_{\rm star}\lt 10^{8.25}\, {\rm M}_{\odot }$|⁠, but as with velocity, it decreases for higher-mass satellites. Although we find little dependence of median d with M_star (not shown), satellites with |$M_{\rm star} \gt 10^8 \, {\rm M}_{\odot }$| today exist at |$\lesssim 250 \, {\rm kpc}$|⁠, whereas there are lower-mass satellites out to |$\approx 800 \, {\rm kpc}$|⁠. Thus, because dynamical friction likely drove higher-mass satellites to smaller velocities and distances, they have smaller orbital lifetimes and smaller ℓ today. Across the full mass range, the mean of the median is |$1.6\times 10^4\, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$|⁠, and the mean range of the 68th percentile is |$0.9\mathrm{ -}2.5\times 10^4\, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$|⁠.

To investigate how ℓ changed over time, we measure the fractional difference in ℓ from first infall into the MW-mass halo to present day, that is, (ℓ(z = 0) − ℓ_infall,MW)/ℓ_infall,MW. Fig. 3 shows this evolution versus |$t^{\rm lb}_{\rm infall,MW}$| (left) and M_star (right). The median ℓ did not significantly change, on average, for satellites that fell in |$\lesssim 10 \, {\rm Gyr}$| ago. However, earlier-infalling satellites systematically increased their angular momenta, by up to a factor of 2 on average, so angular momentum is not conserved for long periods in a dynamic cosmological halo environment, in part because the halo potential evolves on the same time-scale as the orbit. We also stress that the typical 1 − σ width of the distribution is ≈40 per cent, so this represents the typical amount of dynamical scattering that a satellite experienced. Thus, the ensemble satellite population that fell in |$\lesssim 10 \, {\rm Gyr}$| ago did not change in ℓ much, but any individual satellite’s angular momentum changed by up to ≈40 per cent, on average.

Figure 3.

The fractional difference in specific angular momentum since infall into the MW-mass halo, (ℓ_now − ℓ_{infall, MW})/ℓ_{infall, MW}. Lines show the median, and the dark and light shaded regions represent the 68th percentile and full distributions, respectively. Similar to Fig. 3, the vertical grey band shows the free-fall time at the host virial radius. Left: The median fractional difference in ℓ is within ≈10 per cent for satellites that fell in |$t^{\rm lb}_{\rm infall,MW}\lesssim 10 \, {\rm Gyr}$| ago, and the typical scatter about this median is ≈40 per cent. However, earlier-infalling satellites increased ℓ, on average, up to twice their value at infall. Right: Lower-mass satellites show small typical changes in ℓ since infall, though again with a typical scatter of ≈30–50 per cent. At M_star ≳ 10⁸, satellites decreased in ℓ, likely from stronger dynamical friction. Thus, specific angular momentum of satellites is reasonably conserved on average, but with large (≈40 per cent) scatter, for satellites with M_star ≲ 10⁸ that fell in |$\lesssim 10 \, {\rm Gyr}$| ago. Higher-mass satellites, or those that fell in earlier, can experience significant changes in ℓ.

Open in new tabDownload slide

Fig. 3 (right) shows that the median fractional difference in ℓ is minimal for satellites with |$M_{\rm star} \lesssim 10^8\, {\rm M}_{\odot }$|⁠, but again, the 1 − σ width of the distribution is 30–50 per cent. Higher-mass satellites experienced a stronger reduction in ℓ, by roughly 70 per cent, likely from the stronger dynamical friction they experienced.

Fig. 2 thus highlights the strong dependencies of total velocity, E, and ℓ with the lookback time of infall and present-day galactocentric distance, and a lack of dependence on mass, except at |$M_{\rm star}\gtrsim 10^8\, {\rm M}_{\odot }$|⁠. Our results in the middle column especially highlight the large distribution of these properties when selecting satellites at a given distance and, thus, we caution in solely interpreting the median values. Fig. 3 highlights how earlier-infalling and higher-mass satellites experienced larger changes in their specific angular momenta over time, as we explore below.

3.2 Orbital histories

In all but the most ideal conditions, the orbits of satellites change over time. Mechanisms such as the time-dependent (growing) MW-mass host potential, triaxiality of the host halo, dynamical friction, and satellite–satellite interactions can perturb satellite orbits. Because higher-mass satellites experience stronger dynamical friction, and because the MW-mass galaxy and host halo grew over time, a common expectation is that satellite orbits shrink over time, such that their most recent pericentre generally should be their smallest (e.g. Weinberg 1986; Taylor & Babul 2001; Amorisco 2017).

We now explore these expectations and examine trends for the infall times and pericentres in the orbital histories of satellites at z = 0. We investigate ensemble trends to help characterize and compare satellite populations in galaxies such as the MW, M31, and those in the SAGA survey. One should interpret these results for ensembles only, and not necessarily for individual satellites and their orbits, which we will explore further in future work.

Fig. 4 (top) shows the lookback times when satellites fell into the MW-mass halo, |$t^{\rm lb}_{\rm infall,MW}$| (orange) or into any more massive halo, |$t^{\rm lb}_{\rm infall,any}$| (black), as a function of satellite stellar mass, M_star. By ‘any more massive halo’, we specifically mean any halo that is more massive than a given satellite at the same time, which could either be a MW-mass host galaxy’s halo or the halo of a massive central galaxy. We also show a horizontal dotted line to represent when the epoch of reionization ended (e.g. Robertson 2022).

Figure 4.

Lookback times to key events in the orbital histories of satellite galaxies at z = 0 versus their stellar mass, M_star. Lines show the median, and the dark and light shaded regions show the 68th percentile and full distribution across the 13 hosts, respectively. We show a dotted horizontal black line at z = 6 to represent the end of reionization. Top: Lookback time of infall into the MW-mass halo, |$t^{\rm lb}_{\rm infall,MW}$| (orange), and into any (more massive) halo, |$t^{\rm lb}_{\rm infall,any}$| (black). For both infall metrics, higher-mass satellites fell in more recently, with |$t^{\rm lb}_{\rm infall,any}$| ranging from |$\approx 9 \, {\rm Gyr}$| ago for lower-mass satellites to |$\approx 2 \, {\rm Gyr}$| ago for higher-mass satellites, and |$t^{\rm lb}_{\rm infall,MW}$| ranging from ≈7.5 to |$2 \, {\rm Gyr}$| ago. Similarly, the 68th percentile and full distribution for both infall metrics span a similar time range, with the earliest infall times for any surviving satellite reaching nearly |$t^{\rm lb}_{\rm infall,any}\approx 13.2 \, {\rm Gyr}$| and |$t^{\rm lb}_{\rm infall,MW}\approx 12.5 \, {\rm Gyr}$| ago. Satellites with |$M_{\rm star} \lesssim 3 \times 10^{7} \, {\rm M}_{\odot }$| fell into a more massive halo typically |$1\mathrm{ -} 2 \, {\rm Gyr}$| before falling into the MW-mass halo, thus experiencing ‘group pre-processing’. Bottom: lookback time to the most recent pericentre, |$t^{\rm lb}_{\rm peri,rec}$| (green), and to the minimum pericentre, |$t^{\rm lb}_{\rm peri,min}$| (purple), versus M_star, for satellites with at least 1 pericentre. The median in |$t^{\rm lb}_{\rm peri,min}$| decreases from |$4 \mathrm{ -} 5 \, {\rm Gyr}$| ago for lower-mass satellites to |$1\mathrm{ -} 2 \, {\rm Gyr}$| ago for higher-mass satellites. Conversely, the median in |$t^{\rm lb}_{\rm peri,rec}$| is roughly constant at |$1\mathrm{ -}2 \, {\rm Gyr}$| ago across the entire mass range. At most masses, the minimum pericentre occurred |$1 \mathrm{ -} 2 \, {\rm Gyr}$|before the most recent pericentre.

Open in new tabDownload slide

Galaxies form hierarchically, so early-infalling galaxies, either into the MW-mass halo or any more massive halo, were less massive. Additionally, higher-mass satellites experienced stronger dynamical friction, which caused their orbits to lose ℓ and merge with the MW-mass host more quickly. Infall lookback time into any more massive halo decreases from |$t^{\rm lb}_{\rm infall,any}\approx 9 \, {\rm Gyr}$| ago for satellites with |$M_{\rm star} \approx 10^{4.75} \, {\rm M}_{\odot }$| to |$t^{\rm lb}_{\rm infall,any}\approx 2 \, {\rm Gyr}$| ago for |$M_{\rm star} \approx 10^{9.25}\, {\rm M}_{\odot }$|⁠. At |$M_{\rm star} \lesssim 10^7 \, {\rm M}_{\odot }$|⁠, satellites typically fell into another (more massive) halo before they fell into the MW-mass halo. Above this mass, the limited range in mass between the satellite and the MW-mass halo does not leave room for an intermediate-mass host halo. Both the 68th percentile and full distribution for each infall metric span similar ranges; in particular, the 68th percentile ranges from |$4 \,\,\mathrm{ to}\,\, 11 \, {\rm Gyr}$| for satellites at |$M_{\rm star} \lt 10^7 \, {\rm M}_{\odot }$| and |$1 \,\,\mathrm{ to} \,\,9.5 \, {\rm Gyr}$| at higher mass.

The earliest-infalling satellites that survive to z = 0 fell into a more massive halo around |$t^{\rm lb}_{\rm infall,any}\approx 13.2 \, {\rm Gyr}$| ago (z ≈ 8.6) and into the MW-mass halo at |$t^{\rm lb}_{\rm infall,MW}\approx 12.5 \, {\rm Gyr}$| ago (z ≈ 4.7). Furthermore, ≈2/3 of the satellites that fell into their MW-mass halo within the first |$3 \, {\rm Gyr}$| belong to the LG-like paired hosts, presumably because of the earlier assembly of haloes in LG-like environments (see Santistevan et al. 2020). Similar to the analysis of DMO simulations in Wetzel et al. (2015), no surviving satellites at z = 0 were within their MW-mass halo before the end of the epoch of reionization at z ≳ 6 (Robertson 2022). In their analysis, less than 4 per cent of satellites were a member of a more massive halo as early as |$\approx 13.2 \, {\rm Gyr}$| ago (z ≈ 8.6), near to when reionization was ≈50 per cent complete (see e.g. Faucher-Giguère 2020), compared to <1 per cent of the satellites in our sample. Thus, reionization was finished before surviving satellites first became satellites of a more massive halo. For a detailed study on satellite quenching after infall in the FIRE simulations, see Samuel et al. (2022).

Roughly 37 per cent of satellites below |$M_{\rm star}\lt 10^7\, {\rm M}_{\odot }$| were pre-processed before falling into their MW-mass hosts. Using the ELVIS DMO simulation suite of 48 MW-mass hosts, Wetzel et al. (2015) found that for pre-processed subhaloes, the typical halo mass of the more massive halo they fell in was within |$M_{\rm host halo} = 10^{10-12}\, {\rm M}_{\odot }$|⁠, with a median mass of |$10^{11}\, {\rm M}_{\odot }$|⁠. From Fig. 1, this corresponds to a median stellar mass of |$10^{7-9}\, {\rm M}_{\odot }$|⁠. Wetzel et al. (2015) also determined that ∼30 groups contributed all pre-processed satellites, for a typical MW-mass system, with each group contributing between 2 and 5 satellites. In our sample, satellites that were pre-processed typically fell into haloes with |$M_{\rm halo} = 10^{9.3-11}\, {\rm M}_{\odot }$| before falling into the MW-mass host halo, with a median mass at infall of |$M_{\rm halo} = 10^{10.2}\, {\rm M}_{\odot }$|⁠. The stellar masses of galaxies hosted within these more massive haloes ranged from |$M_{\rm star}\sim 10^{6.4-8.9}\, {\rm M}_{\odot }$|⁠, with a median stellar mass of |$M_{\rm star}=10^{7.9}\, {\rm M}_{\odot }$|⁠. Thus, our results are consistent with Wetzel et al. (2015); however, we see a slightly less massive central halo mass.

Fig. 4 (bottom) shows the lookback times to pericentres about the MW-mass host for satellites that have experienced at least one pericentre. Some satellites have orbited their MW-mass host multiple times, so we show the lookback time to two pericentres: the most recent pericentre, |$t^{\rm lb}_{\rm peri,rec}$| (green), and the minimum-distance pericentre, |$t^{\rm lb}_{\rm peri,min}$| (purple). The mass dependence of |$t^{\rm lb}_{\rm peri,rec}$| is weak. However, lower-mass satellites experienced earlier |$t^{\rm lb}_{\rm peri,min}$| than higher-mass satellites, again because lower-mass satellites fell in earlier, when the MW-mass halo was smaller, so their overall orbit and pericentre distance was smaller than for higher-mass satellites that fell in at later times on larger orbits. Given that dynamical friction tends to shrink the orbits of satellites over time, particularly those with |$M_{\rm star} \gtrsim 10^8 \, {\rm M}_{\odot }$| or |${M}_\rm {halo,peak}\gtrsim 3 \times 10^{10} \, {\rm M}_{\odot }$|⁠, one might assume that, for satellites that have experienced multiple pericentres, the minimum pericentre should be equal to the most recent. However, for satellites with |$M_{\rm star} \lesssim 3 \times 10^7 \, {\rm M}_{\odot }$|⁠, the minimum pericentre occurred |$1 \mathrm{ -} 3 \, {\rm Gyr}$| earlier than the most recent pericentre, and the 68th percentile in |$t^{\rm lb}_{\rm peri,min}$| spans a much larger range in lookback time, |$0.25 \mathrm{ -} 9 \, {\rm Gyr}$|⁠, than |$t^{\rm lb}_{\rm peri,rec}$|⁠, |$0 \mathrm{ -} 5 \, {\rm Gyr}$|⁠. At |$M_{\rm star} \gtrsim 3 \times 10^7 \, {\rm M}_{\odot }$|⁠, where the typical number of pericentres experienced is only about 1 (see below), and where dynamical friction is more efficient, the two are more comparable. Thus, the naive expectation that low-mass satellite orbits remain relatively unchanged because they do not experience strong dynamical friction cannot explain the trends in Fig. 4. Furthermore, although the medians between |$t^{\rm lb}_{\rm peri,min}$| and |$t^{\rm lb}_{\rm peri,rec}$| are comparable, the 68th percentile range (and full distribution) shows that differences between these pericentre metrics exist for satellites with |$M_{\rm star} \gtrsim 3 \times 10^7 \, {\rm M}_{\odot }$| also, implying that even the orbits of massive satellites can increase over time. The differences between the most recent and minimum pericentres imply that some mechanism increases the pericentre distances of (especially lower-mass) satellites over time, as we explore below.

Fig. 5 shows trends in |$t^{\rm lb}_{\rm infall,MW}$| (orange) and |$t^{\rm lb}_{\rm infall,any}$| (black) with present-day distance from the host, d. Satellites currently at closer distances typically fell into another halo earlier than more distant satellites. The closest low-mass satellites fell into any more massive halo roughly |$t^{\rm lb}_{\rm infall,any}\approx 10.5 \, {\rm Gyr}$| ago, and this median time since infall decreases to |$4.8 \, {\rm Gyr}$| ago for satellites at d ≳ R_200m. Comparing this to the time since infall into their MW-mass halo, the median |$t^{\rm lb}_{\rm infall,MW}$| is roughly |$0.5 \mathrm{ -} 2 \, {\rm Gyr}$| later across all distances. The range of both the 68th percentile and full distribution generally span similar lookback times, with |$t^{\rm lb}_{\rm infall,any}$| offset to earlier times. The trend of more recent time since infall at larger d is because at earlier times, the MW-mass haloes were smaller and satellites fell in on smaller orbits.

Figure 5.

Similar to Fig. 4 (top), but showing lookback time of infall into the MW-mass halo (orange) and into any (more massive) halo (black) versus present-day galactocentric distance, d. Lines show the median, and the dark and light shaded regions show the 68th percentile and full distribution across the 13 hosts, respectively. Satellites currently closer to the host typically fell in earlier, though with significant scatter. The lookback times of infall for the closest satellites are |$10 \mathrm{ -} 10.75\, {\rm Gyr}$| ago, and these times decrease to |$4.5 \mathrm{ -} 5.5 \, {\rm Gyr}$| ago for the farthest. At any given d today, satellites fell into any more massive halo |$0.5 \mathrm{ -} 2 \, {\rm Gyr}$| before they fell into the MW-mass halo. While d correlates with infall time, we emphasize the large scatter at a given d, which limits the ability of using d to infer infall time for any individual satellite.

Open in new tabDownload slide

Again, one should not interpret our results for individual satellites but rather for populations of satellites. Focusing on the full distribution, which extends across the full range in d, galaxies that fell into their MW-mass hosts between |$t^{\rm lb}_{\rm infall,MW}\approx 3 \mathrm{ -} 9 \, {\rm Gyr}$| ago currently orbit at all distances between |$25 \mathrm{ -} 400 \, {\rm kpc}$|⁠. Thus, although the median shows a clear trend with d, the range in the 68th percentile is |$\gtrsim 2 \, {\rm Gyr}$|⁠, which limits the ability to use the present-day distance of a given satellite to infer its time since infall.

Because of the mass dependence in satellite infall times in Fig. 4, we checked for possible mass dependence of infall time with d by splitting the sample into lower-mass (⁠|$M_{\rm star} \lt 10^7 \, {\rm M}_{\odot }$|⁠) and higher-mass satellites. The difference between |$t^{\rm lb}_{\rm infall,any}$| and |$t^{\rm lb}_{\rm infall,MW}$| exists at all satellite distances in the low-mass sample, and because the stellar mass function is steep, we saw nearly identical results to Fig. 5. However, the higher-mass sample showed little to no difference between the two metrics, with times since infall ranging from |$3.5 \,\,\mathrm{ to}\,\, 7 \, {\rm Gyr}$| ago, because there were not many other more massive haloes for higher-mass satellites to fall into before the MW halo.

Fig. 6 shows trends in the number of pericentric passages, |${N}_\rm {peri}$|⁠, about the MW-mass host (top row) and various pericentric distance metrics (bottom row), versus the time since infall into MW-mass halo, |$t^{\rm lb}_{\rm infall,MW}$| (left), present-day distance from the MW-mass host, d (middle), and satellite stellar mass, M_star (right). Again, we include vertical grey shaded regions that represent the free-fall time at R_200m, t_ff (left column), and MW-mass halo R_200m (middle), as reference values. When presenting trends in |${N}_\rm {peri}$|⁠, we include all satellites, including those that have not yet experienced a pericentric passage, but for trends in pericentre distance we only include satellites that have experienced at least one pericentre. Satellites that have not yet reached first pericentre comprise ≈7 per cent of all satellites.

Figure 6.

For satellite galaxies at z = 0, the number of pericentric passages about their MW-mass host, |${N}_\rm {peri}$| (top), and pericentre distance (bottom) versus their lookback time of infall into the MW-mass halo, |$t^{\rm lb}_{\rm infall,MW}$| (left), present-day galactocentric distance, d (middle), and their stellar mass, M_star (right). We present the median (mean) trends in pericentre distance (number) in the curves, and show the 68th percentile and full distribution via the dark and light shaded regions, respectively. We include all satellites in the top row, but in the bottom row we include only satellites that have completed at least one pericentre. Similar to Fig. 2, the vertical grey shaded bands show the free-fall time, t_ff, at R_200m (left) and the MW-mass halo radii, R_200m (middle). Top row: The mean |${N}_\rm {peri}$| increases with |$t^{\rm lb}_{\rm infall,MW}$|⁠, where the most recently infalling satellites have not had enough time to complete a full orbit (left). This increase of |${N}_\rm {peri}$| for earlier infall is because satellites had more time to orbit, and because those that fell in earlier orbit at smaller distances. Lower-mass satellites experienced more pericentres (right), especially those currently at smaller distances (middle), given their earlier infall and the lack of strong dynamical friction acting on them. Bottom row: Satellites that fell into their MW-mass host earlier experienced smaller minimum pericentres, d_peri,min (purple), and most recent pericentres, d_peri,rec (green; left). Because lower-mass satellites and satellites that are more centrally located fell in earlier, they orbit with smaller d_peri,min and d_peri,rec (middle and right). However, higher-mass satellites also feel stronger dynamical friction that causes them to merge into the host on shorter time-scales. The results highlight how satellite orbits changed between the most recent and minimum pericentres, which often do not occur at the same distance, especially for satellites with |$M_{\rm star} \lesssim 10^7 \, {\rm M}_{\odot }$|⁠.

Open in new tabDownload slide

The top left panel shows the expected trend of more pericentres for earlier |$t^{\rm lb}_{\rm infall,MW}$|⁠. The mean |${N}_\rm {peri}$| is 0 for recently infalling satellites, and it rises to |${N}_\rm {peri}\approx 1$| and is flat across |$t^{\rm lb}_{\rm infall,MW}\approx 2.5 \mathrm{ -} 5.5 \, {\rm Gyr}$| ago, because this time interval is comparable to an orbital time-scale. |${N}_\rm {peri}$| then rises rapidly with |$t^{\rm lb}_{\rm infall,MW}$|⁠, reaching nearly 9 for the earliest-infalling satellites. We compared these trends for lower-mass versus higher-mass satellites and find no significant differences.

Fig. 6 (top middle) shows the dependence of |${N}_\rm {peri}$| on d. Because we find significant differences in |${N}_\rm {peri}$| with M_star (top right), we split the sample into |$M_{\rm star} \lt 10^7 \, {\rm M}_{\odot }$| (solid) and |$M_{\rm star} \gt 10^7 \, {\rm M}_{\odot }$| (dashed). We choose this mass selection given that the lower-mass satellites experience a mean |${N}_\rm {peri}\ge 2$|⁠, while the higher-mass satellites have a mean of |${N}_\rm {peri}\approx 1$|⁠. Lower-mass satellites generally experienced more pericentres than higher-mass satellites at a given d, and the mean number decreases with d for lower-mass satellites from |${N}_\rm {peri}\approx 5$| to 1 for those near R_200m (grey shaded region). Conversely, we do not find dependence on d for higher-mass satellites, with a mean value of |${N}_\rm {peri}\approx 1\mathrm{ -}2$| at all d, likely because of their more recent infall and the increased importance of dynamical friction on them.

Finally, |${N}_\rm {peri}$| declines weakly with M_star, with a mean |${N}_\rm {peri}\approx 2.5$| at |$M_{\rm star} = 10^{4.75} \, {\rm M}_{\odot }$| to |${N}_\rm {peri}\approx 1$| at |$M_{\rm star} \gt 10^9 \, {\rm M}_{\odot }$|⁠. Lower-mass satellites experienced more pericentres, because they fell in earlier (see Fig. 4 top), and also because higher-mass satellites took longer to form and felt stronger dynamical friction that caused them to merge into their MW-mass host on shorter time-scales. Over the full sample, the largest number of pericentres experienced is |${N}_\rm {peri}= 4$| at |$M_{\rm star} \gt 10^7 \, {\rm M}_{\odot }$|⁠, and |${N}_\rm {peri}= 10$| for lower-mass satellites.

The bottom row shows trends for both pericentre metrics: the pericentre with the minimum distance, d_peri,min (purple), and the most recent pericentre, d_peri,rec (green). In the idealized scenario we outline at the beginning of this subsection, an orbiting satellite’s pericentre will remain unchanged or shrink over time because of dynamical friction. However, the panels above show that this often is not true; early-infalling satellites can have larger subsequent pericentres.

The bottom left panel shows the median d_peri,min (purple) and d_peri,rec (green). Both d_peri,min and d_peri,rec are nearly identical for satellites that fell in |$t^{\rm lb}_{\rm infall,MW}\lt 6 \, {\rm Gyr}$| ago, with median values ranging from |$60 \,\,\mathrm{ to}\,\, 100 \, {\rm kpc}$|⁠. For earlier-infalling satellites, both distance trends slowly decrease from |$\approx 100 \, {\rm kpc}$| to |$20 \mathrm{ -} 35 \, {\rm kpc}$|⁠, where the most recent pericentre was roughly |$5 \mathrm{ -} 20 \, {\rm kpc}$| larger than the minimum. Thus, for sufficiently early-infalling satellites that spent longer amounts of time in the evolving MW-mass host halo, the orbits grew slightly over time, which is not expected given the assumption of unchanged orbits or shrinking orbits due to dynamical friction.

We also investigated how the first pericentre a satellite experienced depends on infall time and find qualitatively similar results to the other two pericentre metrics. Earlier-infalling satellites had smaller first pericentres than later-infalling satellites, and the first pericentres were smaller than the most recent ones. The first pericentres were also the minimum a satellite ever experienced in a majority (≈72 per cent) of satellites with |${N}_\rm {peri}\ge 2$|⁠. The only noticeable differences between the first pericentres and d_peri,min occurred for galaxies that fell in |$\gtrsim 6 \, {\rm Gyr}$| ago.

The bottom middle panel shows pericentre distance trends versus current d. As expected, both pericentre metrics increase with d, from |$20 \mathrm{ -} 30 \, {\rm kpc}$| for satellites that are currently closer, to |$70 \mathrm{ -} 80 \, {\rm kpc}$| for satellites near R_200m (grey shaded region). Satellites within |$d \lesssim 225 \, {\rm kpc}$| often had recent pericentres that are larger than d_peri,min by nearly |$10 \mathrm{ -} 20 \, {\rm kpc}$|⁠, so the orbits of these satellites grew. The pericentre metrics in both lower-mass and higher-mass satellites increase with d, but because lower-mass satellites fell in earlier than higher-mass satellites and completed more orbits, they largely drive the differences between d_peri,min and d_peri,rec (and d_peri,min and d_peri,first) in the bottom left panel. We again highlight that the full distributions in both d_peri,min and d_peri,rec span a wide range at a given d, so even though the median trend increases with d, one should not directly apply our results to an individual satellite.

Finally, the bottom right panel shows that lower-mass satellites typically had smaller recent and minimum pericentres. However, at |$M_{\rm star} \gtrsim 10^{8.25} \, {\rm M}_{\odot }$|⁠, the median pericentre distances decrease, likely driven by the onset of efficient dynamical friction. Lower-mass satellites have smaller pericentre distances because they fell in earlier when the MW-mass halo was smaller and less massive. The typical recent/minimum pericentre distance is |$40\mathrm{ -} 60 \, {\rm kpc}$| for satellites with |$M_{\rm star} \lesssim 10^7 \, {\rm M}_{\odot }$|⁠, |$60 \mathrm{ -} 100 \, {\rm kpc}$| for satellites with |$M_{\rm star} \approx 10^{7-8.25} \, {\rm M}_{\odot }$|⁠, and |$\lesssim 100 \, {\rm kpc}$| for higher-mass satellites.

Because the mass of the host can determine the orbits of the satellites, we investigated potential differences between satellites in higher-mass and lower-mass host haloes. At pericentre, satellites are deep in the potential near the galaxy; therefore, the stellar mass of the central galaxy could also correlate with our pericentre metrics. Specifically, we divided the sample in two by selecting the 6 MW-mass hosts with the higher M_star and 7 hosts with lower M_star (see Table 1) and examined their pericentre distances versus |$t^{\rm lb}_{\rm infall,MW}$| and satellite M_star. We find no differences between the two samples versus M_star. Versus |$t^{\rm lb}_{\rm infall,MW}$|⁠, the satellites in higher-mass host haloes had slightly larger, although minimal, d_peri,min and d_peri,rec.

The results in Figs 4–6 suggest a different evolution than expected for some satellites. Lower-mass satellites fell into their MW-mass hosts earlier, when the halo was smaller and less massive, so they complete more orbits than higher-mass satellites in this evolving potential and orbit at smaller distances. Interestingly, the orbits of these lower-mass satellites can increase over time, presumably through the evolving global potential or interactions with other galaxies, which opposes the common expectation of shrinking orbits. However, given the 68th percentile ranges and the full distribution of the pericentre properties, differences between the most recent and minimum pericentres exist at all satellite masses, and not solely at low mass.

3.3 Satellites with growing pericentres

As we showed in Figs 4–6, the most recent pericentre that a satellite experienced is often not the minimum in terms of distance. We now investigate these cases in more detail and refer to satellites with d_peri,min < d_peri,rec as having ‘growing pericentres’.

Satellites with growing pericentres make up 31 per cent of all satellites (ranging from 23 to 46 per cent for a given host). Moreover, growing pericentres comprise the majority (67 per cent across all hosts, ranging from 50 to 86 per cent for a given host) of all satellites with |${N}_\rm {peri}\ge 2$|⁠. In other words, for satellites with two or more pericentres, typically their most recent pericentre was not their closest encounter with their MW-mass host galaxy. Fig. 7 highlights this, showing the fraction of satellites with growing pericentres versus pericentre number. For satellites with |${N}_\rm {peri}\ge 2$|⁠, the growing pericentre population represents >50 per cent of the total sample at any |${N}_\rm {peri}$|⁠, and in some cases, they represent the entire population at a given |${N}_\rm {peri}$|⁠. This fraction broadly increases with |${N}_\rm {peri}$|⁠, at least up to |${N}_\rm {peri}= 6$|⁠, where it represents all satellites, though the fraction fluctuates for |${N}_\rm {peri}$| above that.

Figure 7.

Fraction of satellite galaxies with growing pericentres, relative to all satellites that experienced the given number of pericentres, |${N}_\rm {peri}$|⁠. By definition, these satellites must have |${N}_\rm {peri}\ge 2$|⁠, which is why the fraction is zero for |${N}_\rm {peri}= 0$| and 1. Satellites with growing pericentres represent the majority of all satellites with |${N}_\rm {peri}\ge 2$|.

Open in new tabDownload slide

Although we cannot directly check whether or not this is a temporary occurrence, we compared the most recent pericentre distance to the maximum pericentre a satellite experienced. For satellites that experienced more than three pericentres, we found that roughly 30 per cent of them experienced their maximum pericentre sometime between the minimum and most recent. Of these satellites, the fractional difference between their most recent pericentre distance and their maximum, i.e. (d_peri,recent − d_peri,max)/d_peri,max, was between 1 and 54 per cent, with a median value of 17 per cent, and a 68th percentile range of 8–38 per cent. Thus, because this scenario happens in the minority of satellites, and because the median fractional difference is small, we argue that it is not merely a temporary occurrence. Furthermore, from the top right panel of Fig. 6, satellites with more than three pericentres are generally lower-mass satellites, which will not strongly feel the effects of dynamical friction.

We confirmed that this population of satellites with growing pericentres is not sensitive to the choice for the centre of the MW-mass host in computing satellite distances. Specifically, we examined these trends using the centre of the host DM halo (instead of the centre of the stars in the host galaxy, as is our default). This results in only seven additional satellites whose minimum and most recent pericentres differ by more than 5 per cent, which represents only ≈4 per cent of all satellites with |${N}_\rm {peri}\ge 2$|⁠.

To quantify further the significance of this population, Fig. 8 shows the probability distributions of the difference between key properties at the minimum and most recent pericentres for all satellites with growing pericentres. The left panel shows the fractional difference between the two pericentre distances, (d_{peri, min} − d_{peri, rec})/d_{peri, rec}. As the black point shows, the median fractional difference is −37 per cent, with a 68th percentile range of −15 to −65 per cent.

Figure 8.

Probability distributions of changes in orbital properties for satellite galaxies with growing pericentres. This population comprises the majority (67 per cent) of all satellites that experienced multiple pericentres, or 31 per cent of satellites overall. Each panel shows the median value and 68th percentile via the black square point with error bars. Left: the fractional difference between the minimum and recent pericentre distances, (d_peri,min − d_peri,rec)/d_peri,rec. The median is −37 per cent. Middle: the fractional difference between the specific angular momenta at the minimum and recent pericentres, (ℓ_peri,min − ℓ_peri,recent)/ℓ_peri,recent. The median is −29 per cent. Right: the difference between the lookback times to the minimum and recent pericentre, |$t^{\rm lb}_{\rm peri,min}- t^{\rm lb}_{\rm peri,rec}$|⁠. The median is |$\approx 6.3 \, {\rm Gyr}$|⁠, or 1–2 orbits, given the free-fall time. All metrics show significant differences between the most recent and the minimum pericentre for a given satellite.

Open in new tabDownload slide

Fig. 8 (middle) shows the fractional difference in specific angular momentum at the minimum and most recent pericentres. Nearly all satellites with growing pericentres (>95 per cent) experienced an increase in ℓ between the two pericentres; we do not show in Fig. 8 the small percent with increased ℓ. The median fractional difference in ℓ is −29 per cent, with a range in the 68th percentile of −10 to −60 per cent.

Finally, Fig. 8 (right) shows the difference between the lookback times of the minimum and most recent pericentres. These satellites have a wide range of time differences, with a median of |$\approx 6.3 \, {\rm Gyr}$| and 68th percentile range of |$3.5 \mathrm{ -} 8.5 \, {\rm Gyr}$|⁠. These are slightly longer than the typical orbital periods of these satellites, |$2 \mathrm{ -} 5 \, {\rm Gyr}$|⁠; as Fig. 9 shows, the minimum and most recent pericentres do not always occur successively.

Figure 9.

Orbital distance, d (top), and specific angular momentum, ℓ (bottom), for four representative satellites with growing pericentres, labelled A–D. For each satellite, we list its minimum and recent pericentre distances, along with the fractional difference between the two, and the purple and green arrows show when these events occurred. The top row also shows the growth of the MW-mass R_200m (grey). The majority of satellites with growing pericentres (69 per cent) experienced d_peri,min within |$1 \, {\rm Gyr}$| after infall; in 72 per cent of satellites with growing pericentres, the first pericentre was the minimum, and nearly 71 per cent of this satellite population orbited beyond the MW-mass host R_200m after their first pericentre. The evolution of ℓ shows that the increase in pericentre distance can happen rapidly or gradually, but that it often does not happen near pericentre. Approximately 30 per cent of cases show a sharp increase in ℓ during the first apocentre after infall, which suggests that these satellites may be interacting with other haloes or non-axisymmetric features in the density field.

Open in new tabDownload slide

To provide more context, Fig. 9 shows the orbits (host distance versus time) for four representative satellites with growing pericentres (top row), along with their specific angular momentum (bottom row), labelled A–D from left to right. We chose these four particular satellites at random to span the entire possible range of fractional pericentre distances. The legends show the values of the minimum and most recent pericentres, along with the fractional differences between them; these four satellites range from 12 to 93 per cent. The arrows indicate when these pericentres occurred. For reference, we also show the MW-mass halo’s R_200m(t) (grey).

All four satellites experienced d_peri,min immediately after first infall, and the first pericentre was the minimum in 72 per cent of all satellites with growing pericentres. Furthermore, as Fig. 9 suggests, 71 per cent of satellites with growing pericentres experienced a splashback phase of orbiting beyond R_200m,host after their first pericentre. For comparison, among the population with N_peri > 1 but d_peri,min = d_peri,rec, only 57 per cent experienced a splashback phase. So this suggests that orbiting beyond R_200m,host is associated with a growing pericentre, at least in some cases.

As Fig. 9 (bottom row) suggests, nearly all satellites whose pericentres grew also increased their specific angular momentum. By visually inspecting the histories of the full population, we find that this occurs in two broad ways: (1) steady, gradual increase in ℓ over time, which accounts for 45 per cent of all growing pericentres, and (2) rapid growth in ℓ near a pericentre or apocentre, which accounts for 53 per cent of the satellites. The remaining 2 per cent of satellites are the rare cases in which the pericentres increased from minimum to the most recent, but the angular momentum decreased. The fractional change in ℓ for these satellites is generally small, ≲ 6 per cent. However, some of these satellites show clear signs of interactions with other galaxies, and fell in early ( |$\gtrsim 8.5 \, {\rm Gyr}$| ago). For satellites in category (1), a time-dependent and/or triaxial host halo potential likely plays an important role, especially given that satellites with growing pericentres typically fell in early, when the shape of the host halo potential was changing more rapidly (e.g. Santistevan et al. 2020; Gurvich et al. 2022; Baptista et al., in preparation). Satellite C in Fig. 9 shows a relatively gradual increase in ℓ over time. We defer a more detailed investigation to future work.

For the growing pericentres in category (2), 4/5 of the satellites experienced a rapid increase in ℓ near either a single apocentre, or some combination of them, with the first apocentre being the most common. This is especially apparent in satellites B and D in Fig. 9. The other satellites showed rapid increases in ℓ involving a pericentre that was not the minimum pericentre, much like in satellite A.

Because the fraction of splashback orbits is higher for satellites with growing pericentres compared to the remaining population with multiple pericentres, this suggests that perturbations at d ≳ R_200m,host may play a key role in causing this population. This behaviour is apparent in satellites B and D of Fig. 9, where large spikes in ℓ occur near apocentres, some of which are beyond R_200m,host. These rapid increases in ℓ are caused by rapid increases in the tangential velocities, which typically were of order |$\delta \upsilon \approx 30 \, \rm {km}\, \rm {s}^{-1}$|⁠.

Other satellite mergers with the MW-mass host also can significantly alter the global potential, resulting in orbit perturbations. We investigated correlations of both |$t^{\rm lb}_{\rm peri,min}$| and |$t^{\rm lb}_{\rm peri,rec}$| with the lookback times of mergers, with stellar mass ratios of ≳ 1:100, and did not find a clear correlation between these times. We also investigated correlations of these pericentre metrics with various metrics of host formation times including: the lookback times of when the host galaxy formed 90 per cent of its stellar mass (see Gandhi et al. 2022, for a table of values), the lookback times of when the host formed 10 per cent of its halo mass, and the lookback time of when the host galaxy’s growth transitioned from being dominated by mergers to in-situ formation (Santistevan et al. 2020). We find no significant correlations with these formation metrics.

To investigate whether satellites with growing pericentres have biased orbits, both throughout their history and today, Fig. 10 shows several orbital properties versus stellar mass, for satellites with growing pericentres, all satellites, and all satellites with N_peri > 1 but with shrinking pericentres, i.e. with d_peri,min = d_peri,rec. The top left panel shows the lookback time of infall into the MW-mass halo. Compared to the total sample, as expected, both sub-samples with |${N}_\rm {peri}\gt 1$| fell in typically |$\gtrsim 1 \mathrm{ -} 2 \, {\rm Gyr}$| earlier. However, among the population with |${N}_\rm {peri}\gt 1$|⁠, we find no significant differences between those with growing versus shrinking pericentres, so infall time does not correlate with having a growing pericentre.

Figure 10.

Orbital properties of satellite galaxies at z = 0, including the total sample (blue), those whose minimum and most recent pericentre distances are not the same (‘|${N}_\rm {peri}\gt 1$|⁠, growing pericentres’ in red), and those that experienced multiple pericentres but excluding the ‘growing pericentre’ population (‘|${N}_\rm {peri}\gt 1$|⁠, shrinking pericentres’ in black). Solid lines show median values, and shaded region shows the 68th percentile range for the growing pericentre and total population. Top left: Satellites that experienced multiple pericentres necessarily fell in earlier than the total sample, but among those with |${N}_\rm {peri}\gt 1$|⁠, the growing pericentre population is not systematically biased in infall time. Bottom left: Satellites with growing pericentres experienced |$10 \mathrm{ -} 25 \, {\rm kpc}$| smaller d_peri,min, while the remaining satellites with shrinking pericentres experienced d_peri,min more similar to the total sample. Top right: Because satellites with |${N}_\rm {peri}\gt 1$| fell into the MW-mass halo earlier than the total sample (when the host R_200m was smaller), they have smaller ℓ today, and among these, satellites with growing pericentres have slightly higher ℓ today given their selection. More massive satellites with |${N}_\rm {peri}\gt 1$| have larger ℓ, closer to the total sample, because they fell in more recently. Bottom right: Because satellites with |${N}_\rm {peri}\gt 1$| fell into the MW-mass halo earlier than the total sample and orbit at smaller distances, they are more bound today. Satellites with growing pericentres have orbital energies today comparable to the other satellites with |${N}_\rm {peri}\gt 1$|⁠, though slightly higher at intermediate masses. In summary, satellites with growing pericentres fell in at comparable times to other satellites with |${N}_\rm {peri}\gt 1$|⁠; however, the satellites with growing pericentres orbit with smaller pericentre distances, and with slightly larger ℓ, which makes them unique.

Open in new tabDownload slide

Fig. 10 (bottom left) shows the minimum pericentre distances, d_peri,min. Although all three samples show similar behaviour to Fig. 6, with increasing pericentre distance with increasing M_star, the growing pericentre population is biased to the smallest d_peri,min. The shrinking pericentre sub-sample is generally consistent with the total sample, with typical values spanning |$d_{\rm peri,min}\approx 35 \mathrm{ -} 90 \, {\rm kpc}$|⁠, while d_peri,min for growing pericentre satellites is |$10 \mathrm{ -} 25 \, {\rm kpc}$| smaller, ranging from |$d_{\rm peri,min}\approx 25 \mathrm{ -} 35 \, {\rm kpc}$|⁠. Thus, satellites with growing pericentres orbited closer to the host galaxy. Again, ≈30 per cent of satellites with growing pericentres experienced rapid increases in ℓ during their first apocentre or slightly after d_peri,min. Likely, other important factors contribute to the larger differences between the pericentre metrics, such as the evolving MW-mass host potential, gravitational interactions with other satellite galaxies, and mergers, Fig. 9 hints at, as we plan to explore in future work.

Fig. 10 (top right) shows the specific angular momentum at z = 0. Because the growing and shrinking pericentre sub-samples fell into their MW-mass halo earlier, we expect them to have smaller ℓ. However, the growing pericentres having modestly higher ℓ at z = 0 at most masses, again reflecting that they have scattered to larger ℓ by today.

Fig. 10 (bottom right) shows the specific orbital total energy, E. Consistent with their earlier infall, both sub-samples with |${N}_\rm {peri}\gt 1$| are on more bound orbits today than the total population, at least at |$M_{\rm star} \lt 10^{7.25} \, {\rm M}_{\odot }$|⁠. Any systematic differences between the growing and shrinking pericentres are modest, given the scatter, so we conclude that there are no clear differences in specific orbital total energy at z = 0.

A satellite galaxy can undergo significant mass stripping when it orbits throughout the MW-mass host halo, especially when it is deepest in the host’s potential at pericentre, and this drastic loss in the satellite’s subhalo mass subsequently can affect its orbit. Thus, to better understand the origin of satellites with growing pericentres, including the time-scales over which the orbits changed and potential dynamical perturbations near d_peri,min, we compared the specific angular momentum and DM subhalo mass |$200\, {\rm Myr}$| before and after the minimum pericentre. Near d_peri,min, ℓ changed by a much smaller amount (10−20 per cent) than the change in ℓ from the minimum to most recent pericentres (≈40 per cent). The fractional mass lost near d_peri,min was also minimal (≲ 7 per cent). Thus, in general, the orbital perturbations did not occur just near d_peri,min, as also apparent in Fig. 9.

Finally, we investigated the other orbital properties we presented in the previous figures (total velocity, pericentre lookback times, and the recent pericentre distances) for these sub-samples but find no compelling differences. We find no mass dependence to the fractional differences in pericentre distances or times in Fig. 8, so satellites with growing pericentres exist similarly at a range of masses. Thus, even though higher-mass satellites experience stronger dynamical friction and have smaller orbital lifetimes, we find no mass dependence to whether or not a satellite has an orbit with growing pericentres. We find no strong correlation of the fractional distance or time metrics with either d_peri,min, d_peri,rec, or the lookback times that these occurred. Unsurprisingly, the fractional difference in pericentre distance increased slightly with earlier |$t^{\rm lb}_{\rm infall,MW}$|⁠, given that satellites that orbited for a longer amount of time had more time to experience changes in their orbits.

In summary, of all satellites that experienced |${N}_\rm {peri}\ge 2$|⁠, the majority (67 per cent) experienced a growing pericentre. The most recent pericentre distance is typically ≈37 per cent higher than the minimum experienced, which occurred |$\sim 6 \, {\rm Gyr}$| earlier. Interestingly, about half (45 per cent) of growing pericentres experienced a gradual increase in ℓ, presumably from a time-dependent and/or triaxial MW-mass host potential, and about half (53 per cent) experienced rapid growth in ℓ following either their first or minimum pericentres, during their first apocentres, or during multiple pericentre or apocentre events, which suggests a perturbation by another galaxy. Satellites with growing pericentres are more likely to have been splashback satellites, further suggesting perturbations at large distances. Furthermore, because we measure the orbits of these satellites relative to the MW-mass host galaxy, another effect may be perturbations to the COM of the host galaxy from mergers or massive satellites. Given these complexities and likely multiple causes for the origin of satellite with growing pericentres, we defer a more detailed investigation to future work.

4 SUMMARY & DISCUSSION

4.1 Summary of results

We investigated the orbital dynamics and histories of 473 satellite galaxies with |$M_{\rm star} \gt 3 \times 10^4 \, {\rm M}_{\odot }$| around 13 MW-mass galaxies from the FIRE-2 suite of cosmological simulations. Surprisingly, and in contrast to many (semi)analytical models of satellite evolution, most satellites that experienced multiple orbits experienced an increase in orbital pericentre and specific angular momentum, likely from interactions with the MW-mass host or other satellites. This highlights that satellite orbits do not always shrink and that angular momentum is not always conserved throughout a satellite’s orbital history.

In summary, the topics that we presented in the Introduction and our corresponding results are:

• The relation of orbital properties of satellite galaxies at z = 0to their orbital histories, including lookback times of infall, distances from the MW-mass host, and stellar masses.

  • Satellites that fell in earlier have lower orbital energies and specific angular momenta, though with significant scatter, because satellites that fell in earlier necessarily had to be on smaller orbits to be captured by the MW halo, and the MW-mass host potential continued to grow over time (Fig. 2).

  • Satellites closer to the host generally orbit with higher velocities, smaller specific angular momenta, and have more bound orbits, though with significant scatter (Fig. 2). Total velocity, specific angular momentum, and specific orbital energy do not correlate with M_star except at |$M_{\rm star} \gtrsim 10^8 \, {\rm M}_{\odot }$|⁠, where dynamical friction is more efficient (Fig. 2).

  • Specific angular momentum,ℓ,often is not conserved, even approximately, throughout a satellite’s orbital history (Fig. 3). In particular, earlier-infalling satellites increased in ℓ since infall. More expectedly, higher-mass satellites decrease in ℓ, likely because of dynamical friction. The range of fractional changes in ℓ at smaller M_star and later infall extends ≳ 50 per cent. That said, the average ℓ across the full satellite population remains statistically unchanged since infall.

  • Many lower-mass satellites were pre-processed before becoming a satellite of the MW-mass host. At |$M_{\rm star} \lt 10^7 \, {\rm M}_{\odot }$|⁠, 37 per cent fell into another more massive halo before falling into the MW-mass halo, typically |$\approx 2.7 \, {\rm Gyr}$| before (Figs 4, 5, and A2).

  • No surviving satellites were within the MW-mass halo during the epoch of reionization (z ≳ 6), and less than 4 per cent were satellites of any host halo during this time, similar to Wetzel et al. (2015) (Figs 4 and A2). Surviving satellites at z = 0 fell into the MW-mass halo as early as |$12.5 \, {\rm Gyr}$| ago, and into any host halo as early as |$13.2 \, {\rm Gyr}$| ago.

  • Satellites at a given distance today experienced a large range of infall times into the MW-mass halo. Thus, one cannot infer a precise infall time based solely on a satellite’s present-day distance alone, and the use of total velocity, specific angular momentum, or specific orbital energy alone is similarly limited (Figs 2 and 5).

• Testing a common expectation that the orbits of satellite galaxies shrink over time, that is, that a satellite’s most recent pericentric distance is the minimum that it has experienced.

  • Most satellites at z = 0 with |$M_{\rm star} \lesssim 10^7 \, {\rm M}_{\odot }$| experienced more than one pericentre, while more massive satellites experience only one (Fig. 6), because of their later infall and dynamical friction.

  • Contrary to the expectation that satellite orbits tend to shrink over time, most satellites that experienced two or more pericentres have grown in pericentre distance. Of all satellites with |${N}_\rm {peri}\ge 2$|⁠, 67 per cent experienced a growing pericentre. This represents 31 per cent of all satellites.

  • Typically, the minimum percienter was 37 per cent smaller than the most recent one, because the fractional specific angular momentum increased by 30 per cent (Fig. 8). This minimum pericentre typically occurred |$\sim 6 \, {\rm Gyr}$| before the most recent one (Fig. 8).

  • Satellites with growing pericentres orbited closer to the host (⁠|$d_{\rm peri,min}= 24 \mathrm{ -} 35 \, {\rm kpc}$|⁠) than those with shrinking pericentres.

  • Perturbations at large distances likely contribute to these changes in satellite orbits, given the high fraction (71 per cent) of growing pericentres that were once a splashback satellite. However, we find no single dynamical origin: 53 per cent of satellites with growing pericentres experienced a large increase in ℓ during one or more apocentre, while 45 per cent experienced a gradual, steady increase in ℓ. This suggests that as the MW-mass host halo grows over time, this may help slowly torque the satellites to larger orbits, such that their subsequent pericentres increase. We leave a more detailed investigation of this to future work.

4.2 Inferring infall times from present-day properties

We presented various trends of present-day properties, such as total velocity, specific angular momentum, ℓ, specific energy, E, and distance from the host galaxy, d, with the lookback time of satellite infall, |$t^{\rm lb}_{\rm infall,MW}$|⁠. The median trends in these present-day properties often correlate with |$t^{\rm lb}_{\rm infall,MW}$|⁠. However, we stress that distribution of infall times at fixed property span a large range, limiting the ability to use a property like present-day distance to infer the infall time of a single satellite.

For example, in Fig. 2, while the median specific energy decreases with increasing |$t^{\rm lb}_{\rm infall,MW}$|⁠, for a satellite with a specific energy of |$E = -1 \times 10^4 \, \rm {km}^{2}\, \rm {s}^{-2}$|⁠, the 68 per cent range in |$t^{\rm lb}_{\rm infall,MW}$| is |$1.5 \mathrm{ -} 10.5 \, {\rm Gyr}$| ago. Similarly, although the median specific angular momentum decreases with increasing |$t^{\rm lb}_{\rm infall,MW}$|⁠, a satellite with |$\ell = 2 \times 10^4 \, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$| fell in |$1 \mathrm{ - }9 \, {\rm Gyr}$| ago. Fig. 5 shows that, for a satellite at |$100 \, {\rm kpc}$| today, |$t^{\rm lb}_{\rm infall,MW}\approx 5.5 \mathrm{ -} 10.5 \, {\rm Gyr}$|⁠, and for a satellite near the host virial radius, |$d \approx 300 \, {\rm kpc}$|⁠, it experienced |$t^{\rm lb}_{\rm infall,MW}\approx 2 \mathrm{ -} 8 \, {\rm Gyr}$| ago.

Furthermore, across Figs 2 and 5, at a given satellite total velocity, ℓ, E, or d, the full distribution of infall times spans |$\approx 13 \, {\rm Gyr}$|⁠, nearly the age of the Universe. Thus, while these figures show trends in the median for a population of satellite galaxies, we caution that using any of one of these present-day properties for a single satellite will not precisely determine its infall time into the MW-mass halo. In future work, we will explore how precisely one can infer infall time using full 6D phase-space information, including knowledge about the host potential.

4.3 Comparison to previous work

First, we re-emphasize that these FIRE-2 simulations broadly reflect the observed population of satellites in the LG. Wetzel et al. (2016) and Garrison-Kimmel et al. (2019a) showed that their satellite stellar mass functions and internal velocity dispersions (DM densities) broadly agree with the MW and M31. Samuel et al. (2020) showed that their radial distance distributions broadly agree with the MW and M31, and with MW-mass galaxies from the SAGA survey. Furthermore, Samuel et al. (2021) showed that, although uncommon, spatially or kinematically coherent planes of satellites exist in these simulations, similar to what is observed in the MW and M31. These benchmarks are important for motivating our analysis of their satellite orbits and histories.

Our results agree with Wetzel et al. (2015), who examined similar trends of satellite infall against both M_star and d, using the ELVIS suite of cosmological zoom-in DMO simulations, with abundance matching to assign stellar mass to subhaloes across |$M_{\rm star} = 10^{3-9} \, {\rm M}_{\odot }$|⁠. The mass and spatial resolution in these simulations were |$1.9 \times 10^5\, {\rm M}_{\odot }$| and 140 pc, respectively, ≈5× and ≈3× larger than our baryonic simulations. They found that satellites typically first fell into any more massive halo |$\approx 6.5\mathrm{ -}10 \, {\rm Gyr}$| ago and into the MW-mass halo between |$\approx 5\mathrm{ -}7.5 \, {\rm Gyr}$| ago. These times since infall are consistent with the top panel of Fig. 4 (and the results in Fig. A2) for lower-mass satellites, but the lookback times of infall for satellites at |$M_{\rm star} \gtrsim 10^7 \, {\rm M}_{\odot }$| in our results are more recent than in Wetzel et al. (2015): |$\lesssim 6 \, {\rm Gyr}$| ago for both infall metrics. As we show in Appendix  B, the addition of baryonic physics, especially the additional gravitational potential of the MW-mass galaxy, causes stronger tidal mass stripping and disruption, and because the higher-mass satellites additionally feel stronger dynamical friction, we only see a few higher-mass satellites that happen to survive in our baryonic simulations. As a function of distance from the MW-mass host, Wetzel et al. (2015) also found that satellites experience first infall (into any more massive halo) |$\approx 4 \mathrm{ -} 11 \, {\rm Gyr}$| ago and infall into their MW-mass halo |$\approx 3\mathrm{ -}9 \, {\rm Gyr}$| ago, consistent with our results in Fig. 5.

Rocha et al. (2012) used the Via Lactea II DMO simulation and found a strong correlation between satellite orbital energy and infall time (see also Fillingham et al. 2019; D’Souza & Bell 2022). The authors suggest that satellites that are deeper in the gravitational potential at z = 0 often fell in earlier than satellites farther out and have more negative orbital energies. We find qualitatively consistent values and dependencies of infall time with d as in Rocha et al. (2012): satellites presently closer to the MW-mass galaxy fell in earlier and are on more bound orbits.

Bakels et al. (2021) used one of the N-body Genesis simulations, which have mass and spatial resolution of |$7.8 \times 10^6\, {\rm M}_{\odot }$| and |$1.1 \, {\rm kpc}$|⁠, respectively, which is ≈200× and ≈30× larger than the resolution in our simulations, to study the infall histories of satellites. They analysed the orbits of (sub)haloes of 2309 hosts with |$M_{\rm 200c} \ge 0.67 \times 10^{12} \, {\rm M}_{\odot }$| and found that roughly 22 per cent of all subhaloes are on first infall, much larger than the ≈7–8 per cent in our sample. Furthermore, Bakels et al. (2021) found that roughly 60 per cent of the splashback population of haloes have yet to reach their first apocentre, and the majority (≈86 per cent) have only reached pericentre once, indicating that this population of satellites are on long-period orbits. Given the wide range in MW-mass R_200m, if we select satellites that are currently beyond |$300 \mathrm{ -} 400 \, {\rm kpc}$| to represent the splashback population, we find comparable results: over 95 per cent of the satellites have only experienced one pericentre so far, and the remaining 5 per cent have only experienced two pericentres. For the subhaloes that have experienced pericentre at least once and are currently inside of the host’s virial radius, Bakels et al. ( 2021) found that about half of them have only experienced one pericentre, and ≈30 per cent have not yet reached apocentre. Our result in the bottom left panel of Fig. 6 is generally consistent with this, with a median pericentre number of 1–2 for all satellites in our sample. However, of satellites that experienced at least one pericentre, we find that ≳ 2/3 of them completed more than one, which suggests that the satellites in our simulations completed more orbits. Finally, Bakels et al. ( 2021) noted that roughly 95 per cent of the surviving subhaloes were accreted since z = 1.37 (⁠|$9.1 \, {\rm Gyr}$| ago), where we generally see earlier infall: 95 per cent of our satellites fell into their MW-mass host since z = 2.2 (⁠|$10.7 \, {\rm Gyr}$| ago). Thus, the satellites that survive to z = 0 in our simulations fell in earlier resulting in the larger fraction that have completed more orbits. The differences in the first infall fractions, and the accretion times of satellite galaxies, between our results and Bakels et al. (2021) are likely because of the differences in resolution between the FIRE-2 and Genesis simulations. Because the Genesis simulations have DM particle masses ≈200× larger, they necessarily resolve only more massive satellites.

Fattahi et al. (2020) used the cosmological baryonic zoom-in simulations of MW-mass galaxies from the Auriga project to investigate the |$t^{\rm lb}_{\rm infall,MW}$| for surviving and destroyed low-mass galaxies, and their effect on the growth of the stellar halo. They also found that surviving satellites fell into their MW-mass haloes more recently than the destroyed satellites, similar to the results in Panithanpaisal et al. (2021) and Shipp et al. (2022), who used the same 13 FIRE-2 simulations in our analysis to investigate stellar stream progenitors, their orbits, and their detectability. D’Souza & Bell (2021) also found similar results in their DMO satellite analysis. The analysis by Fattahi et al. (2020) shows similar results to the top panel in our Fig. 4, with more massive satellites falling in more recently. At |$M_{\rm star} = 10^6 \, {\rm M}_{\odot }$| and |$M_{\rm star} = 10^9 \, {\rm M}_{\odot }$|⁠, the authors report average infall lookback times of |$t^{\rm lb}_{\rm infall,MW}= 7.8 \, {\rm Gyr}$| and |$t^{\rm lb}_{\rm infall,MW}= 3.8 \, {\rm Gyr}$|⁠. Our results are broadly consistent, though shifted to more recent times, where satellites at |$M_{\rm star} = 10^6 \, {\rm M}_{\odot }$| and |$M_{\rm star} = 10^9\, {\rm M}_{\odot }$| fell into their MW-mass halo with mean infall lookback times of |$t^{\rm lb}_{\rm infall,MW}\approx 6.5 \, {\rm Gyr}$| and |$t^{\rm lb}_{\rm infall,MW}\approx 2.9 \, {\rm Gyr}$|⁠. We only have three satellites at |$M_{\rm star} \sim 10^9 \, {\rm M}_{\odot }$|⁠, smaller than the sample in Fattahi et al. (2020). The differences between the infall times between satellites in FIRE-2 and Auriga may arise from differences in the stellar mass–halo mass relation at low mass (Grand et al. 2017; Hopkins et al. 2018). Furthermore, the way in which we both average satellite properties over the hosts may contribute to the differences in infall time, given that host galaxies with larger satellite populations will skew the results.

Wetzel et al. (2015) concluded that, in the ELVIS suite of DMO simulations, no present-day satellites were within the MW-mass halo’s virial radius during the epoch of reionization at z ≳ 6. This implies that, for any satellites whose star formation quenched during that time, the MW environment was not the driving factor, so the effects of the MW halo environment and cosmic reionization are separable, in principle. We similarly conclude that no satellites at z = 0 were within their MW-mass halo virial radius during reionization. Although our resolution is still finite, the trend in Fig. 4 is relatively flat with mass, with no indication that it significantly increases for lower-mass satellites. Also, as we show in Appendix  A, we find similar infall trends in subhaloes down to |${M}_\rm {halo,peak}= 10^8 \, {\rm M}_{\odot }$|⁠, which would host ultra-faint galaxies (Fig. A2). Recently, Sand et al. (2022) proposed that the ultra-faint galaxy Tucana B, whose nearest neighbour is |$\approx 500 \, {\rm kpc}$| away, was likely quenched in an isolated environment from reionization. It has an old (⁠|$\approx 13.5 \, {\rm Gyr}$|⁠), metal-poor (⁠|$\rm [Fe/H] \approx -2.5$|⁠), stellar population, and has no recent star formation. Thus, because of its distance to any other massive galaxy, old stellar population, and lack of star formation, Sand et al. (2022) argued that Tucana B is an excellent candidate for a galaxy quenched by reionization. However, our results, and those in Wetzel et al. (2015), imply that no present-day satellites were within a MW-mass halo during reionization. Thus, selecting isolated galaxies today does not necessarily make them cleaner probes of the effects of reionization. Rather, satellites around MW-mass galaxies today provide similarly good candidates to study these effects.

Using the ELVIS DMO simulations, Wetzel et al. (2015) showed that many satellite galaxies first were pre-processed, for |$0.5 \mathrm{ -} 3.5 \, {\rm Gyr}$|⁠, before falling into their MW-mass halo; ≈30 per cent of satellites with |$M_{\rm star} = 10^{3-4} \, {\rm M}_{\odot }$| were members of another group during their infall into a MW-mass halo, and this fraction decreases to ≈10 per cent at |$M_{\rm star} = 10^{8-9} \, {\rm M}_{\odot }$|⁠. Any time before their infall into the MW-mass halo, ≈60 per cent of low-mass satellites were members of another more massive group, falling to ≈30 for high-mass satellites. Over our full sample of satellites, nearly 35 per cent fell into another more massive halo before falling into the MW-mass host, consistent with Wetzel et al. (2015). The fraction of pre-processed satellites in our results is also comparable to the DMO-based results from Li & Helmi (2008) who reported that ≈1/3 of subhaloes were pre-processed, though they selected subhaloes down to |$M_{\rm halo} \gtrsim 3 \times 10^6 \, {\rm M}_{\odot }$| to probe subhaloes that may not host luminous galaxies. Bakels et al. (2021) report that nearly half of all subhaloes with M_sub,acc/M_{host, 200m} ∼ 10⁻³ were pre-processed, and this ratio decreases for increasing subhalo mass. When specifically analysing subhaloes on first infall, Bakels et al. (2021) showed that as many as 40 per cent of subhaloes were pre-processed prior to falling into their MW-mass halo, and this fraction increases for more massive host haloes. More recently, D’Souza & Bell (2021) used the ELVIS DMO simulations to study the times since infall of subhaloes with |$M_{\rm halo,peak} \gt 10^9 \, {\rm M}_{\odot }$| and how they were influenced by a massive merger (>1:10). The distribution of times since infall for their surviving subhaloes range from |$0 \,\,\mathrm{ to} \,\,12 \, {\rm Gyr}$|⁠, and the satellite |$t^{\rm lb}_{\rm infall,MW}$| are peaked towards more recent values compared to the splashback population, which were accreted earlier. The full range of times since infall in our Fig. 4 (and Figs A2 and B1) are consistent with the distribution in D’Souza & Bell (2021). Although D’Souza & Bell (2021) did not specifically focus on the first infall of subhaloes into other more massive satellites/subhaloes, they investigated group infall of satellites and showed that the distribution of time since infall clusters with the timing of the massive merger (and is slightly clustered with lower-mass mergers, >1:15), with many subhaloes becoming satellites of the massive merger |$\lt 2.5 \, {\rm Gyr}$| before it first crossed the MW-mass host radius.

Bakels et al. (2021) showed that after first infall, subhaloes generally lose orbital energy and reach apocentres that are ≈0.8× their turn-around radius, r_ta, and all subsequent apocentres are typically comparable in distance. On the extreme ends, some subhaloes gained or lost orbital energy and, thus, reached larger or smaller subsequent apocentres, respectively, analogous to our satellites with growing pericentres. Regarding the subhaloes that deviate strongly in their first apocentres and r_ta, Bakels et al. (2021) found that nearly ≈2/3 of the satellites with first apocentres ≳ 3r_ta, and ≈80 per cent of the satellites that only reached ≲ 1/4r_ta, were pre-processed. Roughly 1/3 of the satellites with growing pericentres in our sample were pre-processed before falling into the MW, but they may also orbit outside of the more massive halo before falling into the MW halo. Thus, it is unlikely that pre-processing is the only driving factor in the origin and orbital evolution of satellites with growing pericentres.

Both Panithanpaisal et al. (2021) and Shipp et al. (2022) used the same 13 MW-mass galaxies in the FIRE-2 simulations that we use here to investigate stellar stream properties. Stellar streams form via disrupted low-mass galaxies or star clusters; however, before they completely disrupt, because they stretch throughout the halo we learn something about their initial orbits. Shipp et al. (2022) find that systems with smaller pericentres are more likely to form streams, and that the distribution of pericentres in the simulated streams is slightly smaller than the dwarf galaxies in our work. Furthermore, the authors suggest that not only are there differences in the orbital properties of present-day satellites and stellar streams, the orbits of streams with fully or partially disrupted progenitors differ as well, highlighting the complex evolution of low-mass stellar systems.

Finally, D’Souza & Bell (2022) explored uncertainties associated with orbit modelling using the ELVIS DMO simulations. They suggested that using simple parametric models for the MW-mass host (and recently accreted LMC-like galaxy) result in errors that are comparable to the 30 per cent uncertainty in the halo mass of the MW. They also extensively studied the errors associated with modelling the potential of a recently accreted LMC-like galaxy, the initial conditions of the satellites, and the mass evolution of the MW-mass halo, and they show that each comes with errors comparable to or less than the uncertainties in using simple parametric potentials.

Consistent with works like D’Souza & Bell (2022), our results highlight complications and limitations with idealized orbit modelling in a static, non-cosmological MW halo potential; most importantly, our results refute any expectation that the orbits of satellite galaxies always, or even generally, shrink over time. In Santistevan et al. (in preparation), we will use our simulations to pursue orbit modelling of individual satellite histories to compare with idealized orbit modelling in a static host potential.

ACKNOWLEDGEMENTS

We greatly appreciate discussions with Nora Shipp and Pratik Gandhi throughout the development of this paper. We are also thankful for interesting discussion with Andrey Kravtsov on the stellar mass–halo mass relation. IBS received support from NASA, through FINESST grant 80NSSC21K1845. AW received support from: NSF via CAREER award AST-2045928 and grant AST-2107772; NASA ATP grant 80NSSC20K0513; and HST grants GO-14734, AR-15809, GO-15902, and GO-16273 from Space Telescope Science Institute (STScI). JS was supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2102729. RES gratefully acknowledges support from NASA grant 19-ATP19-0068, from the Research Corporation through the Scialog Fellows program on Time Domain Astronomy, from NSF grant AST-2007232, and from HST-AR-15809 from the STScI, which is operated by AURA, Inc., under NASA contract NAS5-26555. We ran simulations using: XSEDE, supported by NSF grant ACI-1548562; Blue Waters, supported by the NSF; Frontera allocations AST21010 and AST20016, supported by the NSF and TACC; Pleiades, via the NASA HEC program through the NAS Division at Ames Research Center.

DATA AVAILABILITY

The python code that we used to analyse these data is available at https://bitbucket.org/isantis/orbit_analysis, which uses the publicly available packages https://bitbucket.org/awetzel/gizmo_analysis, https://bitbucket.org/awetzel/halo_analysis, and https://bitbucket.org/awetzel/utilities. The FIRE-2 simulations are publicly available (Wetzel et al. 2022) at http://flathub.flatironinstitute.org/fire. Additional FIRE simulation data is available at https://fire.northwestern.edu/data. A public version of the gizmo code is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html. Finally, data values in each figure are available at https://ibsantistevan.wixsite.com/mysite/publications.

Footnotes

REFERENCES

APPENDIX A: TRENDS WITH PEAK HALO MASS

In Section 3, we investigated trends of satellite orbital dynamics and histories as a function of satellite M_star, which is the mass most directly observable. However, from Fig. 1, the DM (sub)halo mass of a satellite is 10²–10⁴× larger, so it is the one most important for dynamics. Here we investigate the same trends but as a function of a satellite’s peak halo mass, |${M}_\rm {halo,peak}$|⁠.

We select all subhaloes with |${M}_\rm {halo,peak}\gt 10^8 \, {\rm M}_{\odot }$|⁠, which includes both luminous and dark subhaloes (with no stars). Thus, given the extrapolated abundance matching relations of Moster et al. (2013), Garrison-Kimmel et al. (2017a), and Behroozi et al. (2020) in comparison to our stellar-mass selected sample in Fig. 1, this includes lower-mass subhaloes that likely would host ultra-faint galaxies whose stellar masses our baryonic simulations do not natively resolve. For reference, the fraction of satellites in each mass bin that are luminous in our simulations is: 1 per cent for |${M}_\rm {halo,peak}=10^{8-8.5}\, {\rm M}_{\odot }$|⁠, 14 per cent for |${M}_\rm {halo,peak}=10^{8.5-9}\, {\rm M}_{\odot }$|⁠, 60 per cent for |${M}_\rm {halo,peak}=10^{9-9.5}\, {\rm M}_{\odot }$|⁠, 92 per cent for |${M}_\rm {halo,peak}=10^{9.5-10}\, {\rm M}_{\odot }$|⁠, and 100 per cent above |${M}_\rm {halo,peak}\gt 10^{10}\, {\rm M}_{\odot }$|⁠. Compared to our stellar mass selection, this increases our satellite sample size by a factor of ≈8.5.

Because (sub)halo mass is the most relevant dynamically, but a satellite with a given |${M}_\rm {halo,peak}$| hosts a range of stellar masses given the scatter in the SMHM relation, trends with M_star tend to be noisier. Figs A1, A2, and A3 all show qualitatively similar trends to those in Section 3. In particular, the trends at low halo mass in Fig. A1 show relatively flat dependence, and at |${M}_\rm {halo,peak}\gtrsim 10^{10.5}\, {\rm M}_{\odot }$|⁠, we see a more pronounced decline given the stronger dynamical friction at these masses. Trends with the lookback times of both infall metrics and the pericentre lookback times all qualitatively show similar results and offsets in Fig. A2, and the number of pericentric passages agrees with the stellar mass selection, though it is shifted to slightly smaller |${N}_\rm {peri}\approx 2$| for the smallest subhaloes, compared to |${N}_\rm {peri}\approx 2.5$| at our lowest stellar masses. We do not show trends of pericentre distance, given the lack of a strong dependence on |${M}_\rm {halo,peak}$|⁠, but we compare d_peri,min for satellites in baryonic versus DMO simulations in Appendix  B. In summary, the trends using this halo-mass selected sample are qualitatively similar to the results presented throughout Section 3. Furthermore, our results here imply similar trends for ultra-faint galaxies, where no halo capable of hosting an ultra-faint galaxy, |${M}_\rm {halo,peak}\approx 10^8\, {\rm M}_{\odot }$|⁠, was a satellite of the MW-mass host halo progenitor during the epoch of reionization, z ≳ 6. Similar to the results in our stellar-mass selected sample, we also find that <1 per cent of the satellites in this halo-selected sample were members of a more massive halo during reionization.

Figure A1.

Similar to Fig. 2, the orbital dynamics of satellite galaxies at z = 0, here versus their peak halo mass, |${M}_\rm {halo,peak}$|⁠, for all satellites (luminous and dark). Top: median orbital total velocity is nearly constant with |${M}_\rm {halo,peak}$| at |$\approx 135 \, \rm {km}\, \rm {s}^{-1}$| and 68th percentile ranging from |$71 \,\, \mathrm{ to}\,\, 188 \, \rm {km}\, \rm {s}^{-1}$|⁠. Middle: median specific orbital total energy, E, is nearly constant at |$-0.9 \times 10^4\, \rm {km}^{2}\, \rm {s}^{-2}$| with a 68th percentile range of −2.6 to |$-0.2 \times 10^4\, \rm {km}^{2}\, \rm {s}^{-2}$|⁠. Notably, the dependence on |${M}_\rm {halo,peak}$| is even flatter than with M_star in Fig. 2, though these quantities all decline rapidly at |${M}_\rm {halo,peak}\gtrsim 5 \times 10^{10} \, {\rm M}_{\odot }$|⁠, likely because of the increasing importance of dynamical friction for sufficiently massive satellites. Bottom: median orbital specific angular momentum, ℓ, is nearly constant at |$\ell \approx 1.6 \times 10^4 \, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$| with a 68th percentile range of |$0.9 \mathrm{ -} 2.6 \times 10^4 \, {\rm kpc}\, \rm {km}\, \rm {s}^{-1}$|⁠.

Open in new tabDownload slide

Figure A2.

Similar to Fig. 4, now as a function of peak halo mass, |${M}_\rm {halo,peak}$|⁠, for all satellites (both luminous and dark) at z = 0, showing stronger and smoother trend with |${M}_\rm {halo,peak}$|⁠. Top: similar with the trends with M_star, lower-mass satellites fell into any more massive halo, or their MW-mass halo, earlier than higher-mass satellites, and satellites with |${M}_\rm {halo,peak}\lesssim 10^{10} \, {\rm M}_{\odot }$| typically fell into their MW-mass halo |$\sim 1 \mathrm{ -} 2 \, {\rm Gyr}$| after falling into any more massive halo. As the full distribution shows, no haloes in our sample were satellites of their MW-mass hosts during the epoch of reionization (z ≳ 6), even down to the ultra-faint regime, |${M}_\rm {halo,peak}\approx 10^8\, {\rm M}_{\odot }$|⁠, and less than 1 per cent of satellites were members of another more massive halo during that time. Bottom: As in Fig. 4, the lookback times of both the minimum and most recent pericentres occur earlier for satellites at lower |${M}_\rm {halo,peak}$|⁠. The median |$t^{\rm lb}_{\rm peri,min}$| decreases from |$\approx 4.5 \, {\rm Gyr}$| ago for lower-mass satellites to |$\approx 0.75 \, {\rm Gyr}$| ago at our highest masses. By contrast, the median |$t^{\rm lb}_{\rm peri,rec}$| is much flatter with |${M}_\rm {halo,peak}$|⁠.

Open in new tabDownload slide

Figure A3.

Similar Fig. 6 (top right), now showing the number of pericentric passages about the MW-mass host versus peak halo mass, |${M}_\rm {halo,peak}$|⁠. Now the trend with mass slightly weaker, decreasing from 2 to 1 with increasing |${M}_\rm {halo,peak}$|⁠, because lower-mass satellites fell into the MW-mass halo earlier and on smaller orbits and do not experience as strong of dynamical friction. The full distribution (clipped for visual clarity) is 0–12 pericentres at |${M}_\rm {halo,peak}\approx 10^{8.5} \, {\rm M}_{\odot }$|⁠.

Open in new tabDownload slide

APPENDIX B: BARYONIC VERSUS DARK-MATTER-ONLY SIMULATIONS

Here we compare our results from our FIRE-2 baryonic simulations against satellites in DMO simulations of the same haloes, to understand the effects of baryons and contextualize previous results based on DMO simulations, given that many previous works investigated satellite orbits and infall histories in DMO simulations (e.g. Wetzel et al. 2015; Bakels et al. 2021; D’Souza & Bell 2021; Ogiya et al. 2021; Robles & Bullock 2021), which, among other things, do not model the potential from a central galaxy. Furthermore, stellar feedback in more massive satellites can reduce their inner DM densities, making them more susceptible to tidal disruption (e.g. Bullock & Boylan-Kolchin 2017). With no tidal forces from the central galaxy and with less dense DM cusps within the subhaloes, tidal disruption can be much stronger in baryonic simulations. Recent studies also have used DMO simulations with an embedded disc-like potential (e.g. Kelley et al. 2019; Rodriguez Wimberly et al. 2019; Fillingham et al. 2019; Robles & Bullock 2021).

We compare only simulations that have DMO counterparts at all snapshots, which comprises the 7 MW-mass hosts in isolated environments (names beginning with ‘m12’) in Table 1. As in Appendix  A, for all simulations we select all satellites with |${M}_\rm {halo,peak}\gt 10^8 \, {\rm M}_{\odot }$|⁠, which includes both luminous and dark satellites in the baryonic simulations. In the DMO simulations, we re-normalize |${M}_\rm {halo,peak}$| to account for the loss of baryons by multiplying by 1 − f_b, where f_b = Ω_baryon/Ω_matter is the cosmic baryon fraction. The total number of satellites in the DMO simulations is ≈1.6× higher.

Fig. B1 shows the lookback times of ‘first’ infall into any other more massive halo, |$t^{\rm lb}_{\rm infall,any}$| (top left), specific angular momentum, ℓ (top right), the smallest pericentre experienced, d_{peri, min} (bottom left), and the number of pericentric passages about the MW-mass host, |${N}_\rm {peri}$| (bottom right), for satellites in the baryonic (red) and DMO (black) simulations. Solid lines show the median and the dark and light shaded regions show the 68th percentile and full distribution, respectively.

Figure B1.

Comparing satellite orbital properties versus |${M}_\rm {halo,peak}$| in the baryonic simulations (red) against DMO simulations (black) of the same systems, for satellites at z = 0 within the isolated MW-mass haloes. Top left: satellites in the DMO simulations fell into a more massive halo |$0.5 \mathrm{ -} 3.5 \, {\rm Gyr}$| earlier than in the baryonic simulations. Bottom left: the median minimum pericentre distance, d_peri,min, is smaller in DMO simulations, ranging typically from |$10 \,\,\mathrm{ to}\,\, 55 \, {\rm kpc}$| as compared with |$30 \mathrm{ -} 85 \, {\rm kpc}$| in the baryonic simulations. Top right: the orbital specific angular momentum, ℓ, is smaller in DMO simulations than in baryonic simulations. Bottom right: the mean number of pericentric passages about the MW-mass host is smaller in the baryonic simulations than in DMO. Also, N_peri increases slightly with satellite mass in the DMO simulations, while it decreases monotonically in the baryonic simulations. The primary reason for these differences is that the central galaxy in the baryonic simulations induces stronger tidal stripping and disruption on satellites that orbit near it, leading to a surviving population that fell in more recently, experienced fewer and larger-distance pericentres, and has more orbital angular momentum.

Open in new tabDownload slide

Satellites in the DMO simulations do not feel the gravity of a central galaxy, so they experience weaker tidal stripping and disruption, even if they fell in early or orbit closer to the centre of the halo. Thus, the (surviving) satellites in the DMO simulations generally fell in |$0.5 \mathrm{ -} 3 \, {\rm Gyr}$| earlier than in the baryonic simulations. As a result of the surviving population falling in earlier, satellites in DMO simulations also orbit at smaller distances; they were able to orbit closer to the centre of the MW-mass halo without becoming tidally disrupted, as the bottom left panel shows. Furthermore, surviving satellites have lower ℓ in DMO simulations, given that satellites with smaller ℓ in the baryonic simulations are likely to be tidally disrupted (e.g. Garrison-Kimmel et al. 2017b). Finally, because satellites in DMO simulations fell in earlier and orbit at smaller distances, they completed more pericentric passages (bottom right panel). We also see a small increase in |${N}_\rm {peri}$| with |${M}_\rm {halo,peak}$|⁠, likely because higher-mass satellites in DMO simulations in particular can survive longer than in the presence of a central galaxy.

Our results agree with Garrison-Kimmel et al. (2017b), who compared subhalo populations between DMO and FIRE-2 baryonic simulations using two of the same systems that we analyse (m12i and m12f). They also tested the results of using a DMO simulation with an analytic galaxy potential embedded within the host halo, finding good agreement with the baryonic simulations, which implies that the most important effect in the baryonic simulations is additional gravitational effect of the MW-mass galaxy. They showed that the number of subhaloes between the different types of simulations converges for subhaloes that orbit farther away from the centre of the MW-mass halo. Thus, differences between DMO and baryonic simulations are largest for subhaloes that orbit closer to the centre of the host, where these subhaloes get preferentially disrupted in the baryonic simulations, and result in a satellite population with a larger fraction on more tangential orbits, with higher specific angular momentum.