A promising indirect method to constrain the nature of dark matter (DM) is by measuring the mass function of dark subhalos. At these small scales, different DM models (that have yet to be ruled out observationally) have different predictions for the number of expected subhalos as a function of mass. For example, in the cold DM (CDM) paradigm, the predicted number of dark subhalos with a mass of about 10 7 M e for a galaxy like the Milky Way (MW) is approximately 10 3 subhalos, while in warm DM (WDM) models the expected number depends on the mass of the WDM particles and ranges from 10 to 10 3 dark subhalos (e.g., S. Y. Kim et al. 2018). In self-interacting DM models with velocity-dependent cross section and maximum transfer cross section of about 3.5 cm 2 g -1 , we expect around 100 subhalos (e.g., M. Vogelsberger et al. 2012; V. H. Robles et al. 2019;E. O. Nadler et al. 2021). Thus, measuring the abundance of subhalos at mass scales below 10 7 M e would provide strong constraints between different DM models. A promising way to detect subhalos in the MW is by observing their signatures after they interact with cold substructures such as stellar streams (K. V. Johnston et al. 2002;D. Erkal & V. Belokurov 2015b;J. L. Sanders et al. 2016;J. Bovy et al. 2017;K. Malhan et al. 2021). The prime example is the GD-1 stream, a long (10 kpc, approximately 100°) and thin (20 pc) stream around the MW (C. J. Grillmair & O. Dionatos 2006). The density of the stream is not smooth, but rather shows several gaps in density, one overdensity, and stars orbiting with the stream above the main stream track, known as "the spur" (A. M. Price-Whelan & A. Bonaca 2018; K. Malhan et al. 2019). Detailed orbital modeling has constrained the mass of the perturber to be 10 6 -10 8 M e , possibly originating as a subhalo brought by the Sagittarius dwarf galaxy (A. Bonaca et al. 2019). However, the concentration of this perturber is marginally consistent with predictions for CDM subhalos (T. de Boer et al. 2020;N. Banik et al. 2021bN. Banik et al. , 2021a)). These results have motivated searches for further observational evidence of stream-subhalo interactions across the MW's halo.
Currently, there are about 100 known streams in the MW (C. Mateu 2023) located at different distances and in different regions of the sky, providing a unique opportunity to detect the signatures of dark subhalos across the galaxy. Surveys such as Gaia (Gaia Collaboration et al. 2016), DESI (C. DESI et al. 2016), H3 (C. Conroy et al. 2019), the Vera Rubin Observatory (Ž. Ivezić et al. 2019), WEAVE (G. Dalton et al. 2012), 4most (R. S. De Jong et al. 2019), and Subaru PFS (M. Takada et al. 2014), among others, will further observe streams all the way to the edge of the Galaxy and in external galaxies (S. Pearson et al. 2022;C. Aganze et al. 2023). The data from all of these surveys will provide a multidimensional view (kinematics and chemistry) of the stellar halo that will allow the detection and characterization of the morphology of stellar streams and, hence, the detection of low-mass, completely dark subhalos.
In this section, we detail our approach to studying the stream-subhalo encounter rates in simulations of MW-mass galaxies. We begin by explaining our choice of the simulations and the methodology for identifying the LMC analog. Next, we establish the galactocentric and stream-centric coordinate systems necessary for the analysis. We then outline our subhalo selection criteria and track the LMC-associated subhalos, examining their spatial distribution and assessing their "survival status," i.e., whether they are bound to or disrupted by the LMC analog or the MW.
Furthermore, we introduce our analytical model, outlining the assumptions used to compute the stream-subhalo encounter rates. This involves integrating representative orbits for a suite of both synthetic and real MW streams.
We select two cosmological zoomed-in baryonic simulations of MW-mass galaxies from the Latte suite (A. Wetzel et al. 2023) of the Feedback In Realistic Environments (FIRE) project. 8 These simulations are run with the FIRE-2 physics model (P. F. Hopkins et al. 2018) using the GIZMO code (P. F. Hopkins 2015), which utilizes a TREE+PM solver for gravity and a Lagrangian meshless finite mass solver for hydrodynamics with adaptive spatial resolution. The FIRE-2 model implements star formation and stellar feedback parameters from the Starburst99 stellar evolution models (C. Leitherer et al. 1999) in the ΛCDM cosmology from Planck (P. Collaboration et al. 2016). A detailed description of the FIRE-2 project can be found in P. F. Hopkins et al. (2018), and a description of the Latte suite of MW-mass systems specifically can be found in A. R. Wetzel et al. (2016). D. Horta et al. (2023) made predictions for observable properties of the LMC-analog accretion events in the Latte MW analogs, while N. Panithanpaisal et al. (2021) studied the mass distribution of massive stellar streams and their progenitors. J. Samuel et al. (2020) examined the radial velocity distribution of satellites around isolated and paired MW-M31 analogs and showed that the central disk tidally destroys satellites, altering their radial profile. J. Samuel et al. (2021) showed that spatially thin MW-like planes of satellites can exist in cosmological simulation and the presence of an LMC analog increases the probability of such planes. M. Barry et al. (2023) did a statistical study of spherically averaged subhalo number density and estimated a boost of 1.2-2 due to the LMC analogs. They also estimated that a stream like GD-1 (C. J. Grillmair & O. Dionatos 2006) can have roughly five to six encounters per gigayear.
We select two galaxies labeled m12i and m12b for our analysis. Both are similar to the MW in stellar and gas mass content and their DM mass and density profiles at the present day (P. F. Hopkins et al. 2018;S. Garrison-Kimmel et al. 2018; R. E. Sanderson et al. 2020). Each DM halo has a total mass of about 1. × 10 12 M e , and each simulation uses an initial particle mass of m b = 7100 M e for stars and gas and a DM particle mass m DM = 35,000 M e . m12b includes a massive satellite analogous to the LMC in the MW with its first pericenteric passage about 5 Gyr before the present day, while m12i has a relatively quiescent history with no major mergers for the past 8 Gyr and is used as our "control" system.
We identify the LMC analog in m12b based on the similarity of its orbit around first pericenter (Figure 1) and on the similar merger mass ratio of about 1:8 at pericenter (Table 1). The orbit has d peri = 37.9 kpc and v tot ≈ 350 km s -1 , comparable to the values of d peri = 49 kpc and v tot = 350 km s -1 estimated for the LMC (N. Kallivayalil et al. 2006). Following A. Arora et al. (2022), we define the total mass ratio (TMR) to be the mass of the host galaxy divided by the mass of the satellite at T peri and the pericenter mass ratio (PMR) to be the mass of the MW enclosed within d peri divided by the mass brought in by the satellite within d peri of the center of the satellite. The LMC analog has a TMR of 7.3, while the estimated TMR for the LMC is about ∼4-15, assuming an MW mass of (1-1.5) × 10 12 M e and an LMC mass of (1-2.5) × 10 11 M e (e.g., J. Penarrubia et al. 2016;E. Vasiliev et al. 2021).
The main difference between our analog and the LMC, based on the models in N. Garavito-Camargo et al. (2019), is the plane of the satellite's orbit relative to the host disk and that its first pericentric passage occurred about 5 Gyr before the present day. In our simulation, the LMC analog's orbit is inclined at 107°oriented with respect to the disk plane. However, we do not expect this difference to significantly impact our conclusions. This is because we are interested in computing subhalostream interactions and determining their rates in different halo quadrants relative to the merger trajectory. Moreover, since the orientation of the LMC orbit is not very well constrained owing to uncertainties in measured positions and proper motion and in the MW-mass models (P. J. McMillan 2016; G. Pietrzynski et al. 2019), our simulation with a slightly differently oriented LMC analog provides a plausible scenario. The torque from the disk could affect the orientation and shape of the streams, but the effect is expected to be minor. In fact, a number of studies have investigated the effect of the disk's torque on the shape of tidal streams. For example, D. R. Law et al. (2005) found that the disk torque can significantly affect the orientation and shape of streams for nearly coplanar orbits, while for more inclined orbits the effect is much weaker. On the other hand, S. Kazantzidis et al. (2008) argued that the disk torque is generally not important for the formation of streams, since the subhalo's motion in the host halo is dominated by the overall potential rather than the disk torque. I. B. Santistevan et al. (2024) showed that orbits of satellites in idealized host potentials are completely insensitive to host disk orientation. Furthermore, S. Garrison-Kimmel et al. (2017) showed that tidal disruption of subhalos remains largely unaffected by the geometry of the host disk. Therefore, while the effect of the disk torque on the shape of streams is still debated and likely depends on various factors such as the orbital parameters and the properties of the host and subhalo, it is generally believed to be small for orbits that are not nearly coplanar with the host disk.
The first pericenter passage of the LMC analog in m12b occurs at a look-back time of about 5 Gyr (Table 1), so we define the time T relative to it such that T peri = 0 Gyr. The most massive perturber in m12i reaches the pericenter around the same time and is nearly an order of magnitude less massive relative to its host at the pericenter than the LMC analog we identify in m12b. The age of our host galaxies at pericenter is thus less than the current age of the MW; however, given that the majority of star formation has slowed down by about 7 Gyr in both m12b and m12i (P. F. Hopkins et al. 2018) and the total mass of each halo is close to the MW at the time of the merger, we expect our simulations to be analogous to the real MW even though they are slightly younger than the MW at satellite pericenter. Table 1 summarizes the characteristics of the LMC analog in m12b and the equivalent most massive satellite in m12i at T peri = 0 Gyr. M. Barry et al. (2023) [stop]showed that the presence of an LMC analog can boost the subhalo number density by a factor of 1.2-2. They identified four LMC analogs in simulations m12b, m12f, m12w, and m12c from the Latte suite of FIRE-2 simulations, approximating properties when the analogs approach within 50 kpc of the center and not their actual pericentric distances. This method provides greater statistical accuracy when evaluating the contribution of the LMC analogs to the subhalo population, but not necessarily to their phasespace distribution. Our focus here encompasses a global perspective that also considers the effects of the MW response to the infalling satellites, which is dependent on the satellite's orbit, actual pericenter distance, mass, and velocity. The LMCanalog orbit in m12f is similar to the expected LMC orbit, but the infalling satellite is only half as massive, with a TMR of about 16 (A. Arora et al. 2022). Consequently, the host's response to the LMC analog in m12f is significantly weaker.
Moving on to m12w and m12c, both simulations feature satellites similar in mass to the LMC. However, m12w stands out owing to its highly eccentric, fast-moving orbit, resulting in a very radial merger with its first pericenter at 8 kpc and eventual complete tidal disruption of the satellite. It is almost 6.5 Gyr before the present day (1.5 Gyr earlier in comparison to m12b). This makes it highly unreliable for predictions. In contrast, the LMC analog in m12c is closest to the present day but follows an orbit completely within the disk plane, with the actual first pericenter occurring at 18 kpc. This particular orbit is not expected to induce a north-south asymmetry. Selecting an arbitrary pericenter at 50 kpc would not accurately represent halo deformations. For the sake of completeness, we have included our encounter rate calculation for m12c in Appendix C.
The MW-LMC system we consider involves two galaxies orbiting around their common COM, which is moving with relatively constant velocity through space in a cosmological simulation box. The frame centered on and moving with the system COM is thus an inertial frame (modulo interactions with the next most massive satellite galaxy of the MW). We will refer to this as the system frame, marked by primed quantities, related to the one in which the simulation is run (the "simulation box frame").
From our perspective in the MW's disk, however, the Sun orbits the center of the MW (the local minimum in the potential and maximum in density), and we commonly use a coordinate system centered on and moving with this location, referred to as the galactocentric frame. This frame is not inertial, since the MW and the LMC are orbiting their common COM, but instead is related to the system frame by
where the unprimed quantities are in the galactocentric frame and x MW ¢ and v MW ¢ denote the position and velocity of the MW's center in the system frame, respectively, both of which are functions of time (as explicitly highlighted in this equation).
Since v d dt 0 MW ¹ , the galactocentric frame is not inertial. The displacement of the galactocentric frame in position and velocity is defined so that at the time of pericenter t peri (the present configuration in the actual MW-LMC system)
where x ref ¢ is the distance between the MW center and the system COM and v ref ¢ is the relative motion of the MW's center induced by the LMC, also known as the "'reflex motion," which for the real MW-LMC system is approximately 30-40 km s -1 at the location of the Sun (D. Erkal et al. 2019; M. S. Petersen & J. Peñarrubia 2020; E. Vasiliev et al. 2021;E. Vasiliev 2023).
The LMC also accelerates the subhalos orbiting in the MW globally, leading to a net bulk motion of these subhalos directed toward the MW-LMC center. The acceleration experienced by an individual subhalo varies depending on its location in the galaxy with respect to the LMC. The LMC would roughly equally accelerate both the subhalos around a stream and the stream itself. However, if subhalos are anisotropically or asymmetrically distributed around the stream, small variations in the accelerations experienced by the subhalos due to the LMC will lead to a perturbatively induced nonzero bulk motion of the subhalos around the stream. In contrast, if the subhalos are isotropically distributed, the net bulk motion around a stream will be zero.
Given that the simulations are run in an arbitrary simulation box frame, the traditional approach is to establish a "principal axis" by aligning the galactic disk to the X-Y plane at the time of the analog's first pericentric passage. In this frame, the LMC analog in m12b is located at ( ) x 23.3, 26.6, 13.7 analog = -kpc at T peri . However, we establish an additional "rotated axis," which aligns the position unit vector of the LMC analog in m12b at its first pericentric passage with the position unit vector of the real LMC in the MW at its first pericentric passage, located at x LMC = ( ) 2.3, 20.2, 41.1 -kpc, based on the orbits presented in N. Garavito-Camargo et al. (2019). Details on how to compute this rotation can be found in Appendix A.
Figure 1 plots the trajectory of the LMC analog (red/green) in m12b in the rotated axis galactocentric coordinates, along with the total distance from the MW's central location (on the bottom right), compared to an MW-LMC constrained simulation (blue) from N. Garavito-Camargo et al. (2019). The green sections of the trajectories indicate the periods over which we compute encounter rates. The overall orbit and the present-day location of the analog match reasonably well those of the LMC. Unless otherwise noted, all analyses and results in this paper are presented in the rotated axis frame.
We will compute the subhalo-stream encounter rates for streams distributed in four quadrants in the sky. We define each quadrant in the galactocentric-rotated axis frame in Table 2 for reference, also marked in the Figure 1 Y-Z panel.
Finally, to compute quantities relevant to stream-subhalo encounters, the relevant frame is the one moving with the stream, here referred to as the stream-centric frame. Since most of the currently known streams are closer to the MW center than the LMC, we will relate the stream-centric frame, denoted by x and ṽ, to the galactocentric frame x v , (for streams orbiting beyond the LMC one would substitute the system frame x¢): Note. T peri : time of closest approach ("pericenter") between the main galaxy and the satellite (T = 0 Gyr). All properties are evaluated at T peri . M å,host : stellar mass of the halo. d peri : pericenter distance between the satellite and the main galaxy. TMR: total mass ratio of the MW and the satellite, M main /M sat . PMR: pericenter mass ratio, M main (<d peri )/M sat (<d peri from sat ). v rad , v tan : radial and tangential velocities of the satellite with respect to the MW. Parameter v rad is nonzero owing to finite time resolution between snapshots around the pericentric passage.
where
are the position and velocity, respectively, of a representative particle in the stream in the galactocentric coordinates (in other words, the stream's representative orbit). The relevant quantity for determining the rate of subhalo interactions is the distribution of the subhalos velocities perpendicular to the stream, which corresponds to the cylindrical radial velocities of subhalos in this stream-centric frame ṽR .
In order to compute the encounter rates, we need to inject stellar streams into the simulations and integrate their orbits; here we describe the methodology we employ.
Each stream is represented as a test particle that is integrated using the time-dependent potential of the system. We generate two sets of stream catalogs: one resembling the MW streams, and one a synthetic catalog of 5000 streams. The properties and the initial conditions (ICs) used to generate these two catalogs are described below:
1. MW streams: 97 streams orbiting the MW with phasespace properties resembling the observed streams of the MW. We use galstreams (C. Mateu 2023) to compute the median phase-space values of each stream. We then inject the stream's positions and velocity into our rotated axis frame, aligning the orientation of the streams and the LMC analog in our model with that of the real LMC. 2. Synthetic streams: We generate 5000 synthetic stream test particles uniformly spread across the sky in a Mollweide projection and different orbital properties.
We simulate these streams to be at their pericenter, with distances between 20 and 80 kpc from the galactic center at T = 0 Gyr and varying azimuthal velocities and inclination velocities such that the total tangential velocities at T = 0 Gyr are between 250 and 400 km s -1 (set by the real MW streams). These initial conditions produce a uniform distribution of streams' orbital eccentricities and pericentric distance in all quadrants.
These ICs are used to start the orbital integration for each stream at T = 0 Gyr. We then integrate the orbits for a total of 0.5 Gyr, corresponding to the shaded green band in Figure 1 around the pericenter of the LMC-like satellite. We integrate backward for -0.4 Gyr starting at T = 0 Gyr and then forward by 0.1 Gyr starting from T = 0 Gyr. A value of 0.5 Gyr is sufficient to induce noticeable morphological changes in a stream's structure, such as the emergence of open gaps and the occurrence of kinks (D. Erkal & V. Belokurov 2015a). In the following section, we describe the time-dependent model that we use in our orbital integration.
To integrate the orbits of the test particles (each representing a stream), we use time-evolving potential models from A. Arora et al. (2022) following the scheme described in A. Arora et al. (2024). In short, we fitted the potentials with basis function expansions (BFEs) on the host halo at every snapshot in the simulations within the 20 snapshots that span the 0.5 Gyr around the T = 0. These models can adequately describe deformations in the halo caused by the LMC analogs (e.g., M. S. Petersen & J. Peñarrubia 2020;N. Garavito-Camargo et al. 2021;A. Arora et al. 2022; M. S. Petersen et al. 2022). BFEs have also been shown to successfully reproduce orbits, even in the presence of massive satellites, for short periods of time (less than 1 Gyr) in both idealized simulations (E. Vasiliev 2024; S. Lilleengen et al. 2023) and cosmological simulations (B. Lowing et al. 2011;J. L. Sanders et al. 2020;A. Arora et al. 2024).
Since we represent each stream by the orbit of a single test particle, we do not simulate the perturbations along the leading and trailing arms of the stream (as is observed, e.g., in the reallife Orphan stream; S. Lilleengen et al. 2023).
The integration is performed in the principal axis frame on which the BFE models are fit. Nonetheless, we report our results in the rotated axis frame. Our orbit integration of subhalos and streams is not entirely self-consistent, as the subhalos respond to the LMC analog in a live simulation, while the streams are integrated within a smooth, time-evolving potential model. However, one can simply integrate the subhalo orbits within the time-dependent potential using a prescription to account for tidally destroyed subhalos and count flybys. This approach is certainly feasible, but it also relies on the assumption that subhalo orbits can be properly reproduced, without taking into account the effects of dynamical friction.
We use the integrated orbits to calculate the apocenter and pericenter distances, and we categorize the streams as prograde or retrograde relative to the orbit of the LMC analog in m12b. Only about 13 of the known MW streams have pericenter distances greater than 20 kpc from the galactic center in the m12b potential.
We use the ROCKSTAR halo finder (P. S. Behroozi et al. 2012a) to identify DM subhalos. Merger trees were constructed to link the catalogs at each snapshot using the consistent-trees code (P. S. Behroozi et al. 2012b), as detailed in J. Samuel et al. (2021) and N. Panithanpaisal et al. (2021). Using these halo catalogs, we identify subhalos that are within the MW and brought in by the LMC analog. We impose a mass range of 10 6 M e M sub 10 9 M e on all the subhalos based on our resolution. We select the subhalos bound to and within the virial radius of the LMC analog 2 Gyr before the first pericenter, when the d LMC + R vir,LMC is greater than the R vir,MW , and track their evolution through time until disruption using consistent-trees.
The LMC analog brings a total of 1029 subhalos. Figure 2 plots the location of the LMC-analog subhalos as density contours in the X-Z plane for m12b at -1 Gyr (left), 0 Gyr (middle), and 1 Gyr (right) with the LMC-analog trajectory and We track whether these subhalos are bound to the LMC analog or the MW using a total negative energy criterion (i.e., E sub,tot 0 in the LMC analog or the MW reference frame). We also note whether the subhalos merge into the LMC analog or the MW. A subhalo is considered "merged" when consistent-trees can no longer track its center and its descendant subhalo links to either the MW or the LMC analog. consistent-trees can successfully track subhalos with masses greater than 10 6 M e (J. Samuel et al. 2021).
Figure 3 illustrates the number of LMC-analog subhalos bound to the LMC analog (blue) or the MW (orange), merged into the LMC analog (green) or the MW (red) as a function of time starting from 2 Gyr before the first pericenter passage. By the first pericenter (0 Gyr), 80% of the subhalos are tidally stripped from the satellite and bound to the MW. By 1 Gyr, only half of the subhalos survive the merger, while the rest are destroyed and form dark streams. In the scope of our simulations, and given our inability to resolve these dark streams, we leave them out of our calculations for future research. These dark streams are less likely to leave observable signatures after their interactions with stellar streams because of their low surface density.
Counting the number of stream-subhalo encounters is not feasible, due to the limited temporal resolution in the simulations. Our time range consists of approximately 20 snapshots, each spaced 25 Myr apart. Therefore, we calculate the encounter rate at each snapshot from our simulation between the subhalos and a stellar stream using the counting scheme described by D. Erkal et al. (2016), which is similar to As the analog falls in the MW, the LMC-analog subhalos are tidally stripped away from it and settle into bound orbits around the MW. This is obvious at 0 and 1 Gyr, as the density contours are more widespread in both cases. At 0 Gyr, the majority of the subhalos lag behind the infalling LMC analog, while some completely settled into newer orbits, producing a dipole term in density. At 1 Gyr, almost all of the subhalos are dispersed away from the LMC analog, as the central density of subhalos is away from the LMC analog. Also note that there is a slingshot action as a few subhalos end up unbound and outside the virial radius of the MW after the analog's flyby (upper right corner for T = 1 Gyr). the one conducted by J. H. Yoon et al. (2011). The encounter rates are computed as the number of subhalos entering a cylinder of radius b max (impact parameter) around a straight line stream of length ℓ s . The encounter rates depend on the local number density of subhalos n sub around the stream and the cylindrical radial velocity ṽR distribution of subhalos with respect to the stream.
The number of subhalos passing through a cylinder of length ℓ s and radius b max at time t snap , within a time interval dt, is then given by
Here ( ˜) P v R represents the probability distribution function of ṽR . In contrast to D. Erkal et al. (2016), where they considered a straight line stream of length ℓ s and computed encounter rates by modeling a cylinder of length ℓ s and radius b max around the stream, our approach simplifies the stream to a single-particle representation. This approximation does not accurately describe a cylinder around the stream. Instead, we calculate ṽR as the relative velocity component perpendicular to the tangent plane of the stream at xstr as ˜(
, with direction toward the plane, from both the top and bottom of the plane, represented with negative values. This approach overlooks the orientation of the stream on the tangent plane itself. Subhalos moving toward the plane with a zero component perpendicular to the plane can still perturb the stream, an effect not accounted for in our single-particle model. Therefore, we integrate over different orientations of the plane assuming that ( ˜) P v R remains symmetric under small rotations, imitating the change in direction of the normal to the cylinder. This produces an extra factor of 2π. In this calculation, we consider only subhalos entering through the sides of the "cylinder," so we integrate over all negative ṽR . We can then rewrite the encounter rates as
where ( ˜) v R is the first moment of the probability distribution of subhalos entering the cylinder ( ṽ 0 R ). ( ˜) P v R is often approximated by a Gaussian distribution of ṽR assuming a nonzero mean ( m) and dispersion (s) (J. H. Yoon et al. 2011;D. Erkal et al. 2016). In such a scenario, ( ˜)
. Parameter ( ˜) f g is a unitless scaling factor that depends on the mean and dispersion of ṽR , hereafter the anisotropic boost factor. This factor equals unity when ˜0 m = (see Figure 15 in Appendix B for more details) for the traditional encounter rate formulation such as in J. H. Yoon et al. (2011), D. Erkal et al. (2016), and M. Barry et al. (2023).
In Appendix B, we demonstrate that this approximation is valid to within 10% error compared to the numerical integration.
The density of subhalos, n sub , is computed as a function of time by counting the number of subhalos within a 10 kpc width spherical shell centered around the stream. This shell is further confined to a sky slice that spans ±π/2 in colongitude and ±π/4 in colatitude around the stream's position. The stream is always positioned at the center of this slice at each time step. We then divide the counts by the volume of this sliced shell, which is one-fourth of the total volume of the spherical shell. Our decision not to use a local volume around the stream was made to mitigate numerical noise owing to low subhalo counts. Similarly, we evaluate ( ˜) P v R from subhalos within the same shell and in the stream's quadrant for each time step of the stream. We numerically integrate Equation (9) using ( ˜) P v R to compute the encounter rates for each time snapshot using Equation (8).
Finally, we sum the number of encounters per gigayear for a stream using the equation
where Δt snap represents the time interval between snapshots, the sums are taken over all time steps t snap , and dN dt enc is computed using Equation (8).
In this paper we do not explore the impacts caused by any single encounters in our analysis, which are more sensitive to the individual stream and subhalo kinematics and require a case-bycase study. We rather focus on a global analysis of the encounter rates in this paper. Traditionally an impact weighting based on the impulse approximation
~, where b is the distance of the closest approach between a subhalo and a stream, has been used in the literature (D. Erkal & V. Belokurov 2015a; D. Erkal et al. 2016). However, implementing such weighting in our method would require integrating over the mass distribution of subhalos, which is trivial since P(M sub ) ∝ 1/M sub , and a distribution of b. Finding b requires computing pairwise distance between each stream and all the subhalos within some set distance cutoff around the stream during the integration time, and the number of subhalos changes nontrivially. Alternatively, one can assume that b is a nonlinear function of the velocity perpendicular ( ṽR ) and parallel ( ˜|| v ) to the stream. These assumptions lead to a coupled 2D integral with both the velocity distributions in the denominator and ˜˜|| v v R 2 2 + in the numerator; this integral diverges at 0 and is highly sensitive to minor fluctuations. However, in our case the impact weighting for the subhalos brought in by the LMC analog will be similar to the impact weighting of the MW subhalos; this is demonstrated by showing that the velocity field is fairly independent of the subhalo mass in Section 3.
In this section, we illustrate how and why each of the physical quantities in the encounter rates (subhalo number density and radial velocity distribution) is affected during the infall of the LMC analog. In Section 3.1 we start by exploring the temporal evolution of the phase-space distribution of the host and contributed subhalos. In Section 3.2 we explore the temporal changes in the radial velocity distribution of the subhalos caused by the perturbations of the LMC analog. In Section 3.3 we explore the radial and angular variations on the sky in the integral of the probability distribution ( ) v rad . Lastly, we quantify the evolution of the number density of subhalos in Section 3.4.
We begin by characterizing the phase-space distribution of both the host subhalos (those within the virial radius of the MW and not brought in by the LMC analog) and the subhalos brought in by the LMC analog, hereafter the contributed subhalos. We show how the presence of these subhalos influences the overall phase-space distribution and dynamics of the MW subhalos.
Figure 4 shows the phase-space distribution at three equally spaced time steps: 1 Gyr before pericenter (left column), when the LMC analog is at pericenter (middle column), and 1 Gyr after pericenter (right column). The ratio of contributed subhalos to the total subhalos as a function of radius is plotted in the top row. Despite the contributed subhalos composing only a fraction (less than 10%) of the host subhalos, they can increase the spherically averaged subhalo density near the LMC analog's radial distance (indicated by the red dashed line) by up to 20%.
The middle and bottom rows show the total velocity and radial velocity distributions in the galactocentric coordinates of the contributed subhalos (depicted as open contours), relative to the background MW subhalo distribution (depicted as filled contours) within the virial radius of the MW. Notably, both the distributions of the contributed subhalos within the virial radius reveal a distinct leading arm (higher density closer to the host center) compared to the LMC analog (indicated by the red star) and exhibit different phase-space characteristics compared to the MW. As the merger progresses, the surviving subhalos gradually phase-mix into new orbits resembling those of the MW subhalos, which becomes evident as early as 1 Gyr after the first pericentric passage, aligning with findings in L. V. Sales et al. (2017) and M. Barry et al. (2023). The tangential velocity distribution of the subhalos (not shown in the figure) in the galactocentric coordinates mirrors the trends seen in the total velocity distribution. The escape velocity curves are overlaid in the middle row at each time step. Most of the subhalos are below the escape velocity; however, some exhibit higher total velocities. Many of these subhalos are falling toward the galactic center and thus experience tidal stripping. The insets in the upper right corner of the middle row show the mass distribution of the contributed subhalos within the virial radius.
It is important to highlight that the LMC analog's position is not centered within the open contour distribution. This is because we are solely analyzing the subhalos situated within the MW's virial radius.
Furthermore, at T = -1 Gyr, the mean radial velocity of the contributed subhalos is ≈150 km s -1 (marked by the green dashed line), corresponding to the analog's radial velocity. By T = 0 Gyr, the mean radial velocity of the surviving contributed subhalos aligns within ±25 km s -1 of the MW subhalo population as they experience phase mixing.
In summary, the presence of contributed subhalos can increase the spherically averaged subhalo density by up to 10%-40% (Appendix C); M. Barry et al. (2023) note a higher boost due to the contributed subhalos in other LMC analogs. This effect will become more pronounced when considering the azimuthal dependence of the subhalo density distribution. Surviving LMC-analog subhalos gradually phase-mix into orbits resembling the MW subhalos. However, before T = 0 Gyr, the contributed subhalos can shift the mean radial velocity of subhalos, thus affecting the stream-subhalo encounter rates (Equation ( 8)). These effects are most prominent during the first pericentric passage and can either increase or decrease the encounter rates.
The local subhalo density and velocity distributions are not only affected by the contribution of additional subhalos from the LMC analog; they are also perturbed by the shift of the COM frame induced by the presence of the LMC analog. To illustrate this, we show in Figure 5(a) the COM velocity of the MW's host over time for both m12b (blue) and m12i (green) in the reference frame of the cosmological volume. The MW COM is determined using the shrinking spheres method described in C. Power et al. (2003). The total central velocity is calculated using the mean velocities of stars within 10 kpc of the MW's COM.
As the LMC analog approaches its first pericentric passage, there is a 60 km s -1 change in the MW central velocity primarily caused by the LMC analog moving the COM of the system away from the host's center. This persists for 1-1.5 Gyr after the first pericenter as the LMC analog and the MW orbit their common COM. The effect reaches another peak as the system approaches its second pericenter passage at approximately 1.7 Gyr. In contrast, such motion is not observed in m12i, which has a relatively constant velocity of about 130 km s -1 in time with deviations of order 10 km s -1 , consistent with the results from J.-B. Salomon et al. (2023).
In As the LMC analog transitions from the Q.III to the Q.IV and moves closer to the galactic center, the MW COM moves toward the trajectory of the LMC analog. This effect is most pronounced when the LMC analog approaches its first pericenter at T = 0 Gyr. By the time of the LMC analog's second pericenter, the induced reflex motion is not as strong, as the halo starts to relax and the satellite loses energy. In contrast, m12i does not show a systematic trend with the merging satellite and exhibits a random distribution of velocity vectors. The green cross represents the average central displacement over the duration of the satellite merger.
We now explore how the velocity of the COM of the host changes the radial velocity of the subhalos as a function of time. In Figure 6, we plot the mean radial velocity of subhalos in galactocentric coordinates 〈v rad 〉 in multiple 20 kpc radial bins. We observe variations in the mean radial velocities at the LMC analog's first pericenter (T = 0) and in subsequent apocenters and pericenter in m12b for all the subhalos within the virial radius (solid black). These variations are due to both the motion of the COM of the host and the contribution from the LMC-analog subhalo population with its own distinct phase-space properties. Similarly, the mean tangential velocity of subhalos in the galactocentric coordinates also changes by about 30 km s -1 . The most significant deviations are observed in the 40-60 kpc bin (solid green) and the 60-80 kpc bin (solid blue), corresponding to the distances probed by the LMC analog during the given time interval. The high variations in the 20-40 kpc bin are attributed to a small number of statistics. m12i (dashed line) demonstrates a relatively low and constant mean value. Interestingly, the velocity dispersion (not shown here) remains relatively constant azimuthally at a radius and in time in both m12i and m12b.
The changes in the velocity distribution of the subhalos also affect the first moment of the velocity distribution ( ) v rad in the galactocentric coordinates. Here we quantify those changes as a function of time. Figure 7 shows the temporal evolution of ( ( ) v rad ) computed numerically using Equation (9) for subhalos in different mass bins (color-coded). These subhalos orbit within a distance range of 40-80 kpc from the galactic center. The ( ) v rad for subhalos in different mass ranges exhibits similar trends over time (varying between 50 and 100 km s -1 ), yet the amplitudes are different by 25%. Notably, there is an increase in ( ) v rad right before the LMC analog's first pericentric passage at 0 Gyr as a result of the COM motion of the host. The black curve, which accounts for all subhalos, closely matches the lower mass bin curves, reinforcing that the overall trend is consistent. The high jitter observed in the 10 7.5 -10 9 M e mass range results from small number statistics, which causes fluctuations in the integral and affects its convergence owing to the lower number of subhalos in this mass bin. The consistency in ( ) v rad highlights that the overall velocity field is mostly decoupled to the mass and density field of subhalos. While we have only considered spherically averaged effects here, this decoupling between the velocity field and the mass of the subhalos is expected to persist even in spatially localized regions.
As we have demonstrated in Section 3.1, the location of the contributed subhalos from the LMC analog varies as the satellite orbits the host. Moreover, due to the host's COM motion, the radial velocities of the entire population of subhalos change as a function of radius and time. Here we further explore the angular variations that the radial velocity distribution induces in the ( ) v rad . To simplify our calculation, we assume that the radial velocity distribution around the streams is Gaussian. We found this approximation to be accurate at the 10% level as shown in Appendix B. However, we do not assume that the mean velocity of the radial velocity distribution m is zero. To quantify the degree of asymmetry in different regions of the sky, we compare our results to the case of ˜0 m = .
In Figure 8, the top panel plots the LMC-analog distance from the center as a function of time, with the LMC analog's quadrant marked with shaded color regions at each time step. We also show selected 20 kpc bins within the 20 and 80 kpc interval, based on the LMC analog's location. In the middle panel, we plot the change in the anisotropic boost factor, that is, [f (γ)f (0)], as a function of time in the selected 20 kpc bins (show in the top panel) for m12b with contributed subhalos (labeled as m12b; dark blue) and without contributed subhalos (labeled as m12b (without LMC subhalos); magenta) and m12i (green), the galaxy without massive satellites. The bottom panel shows the change in the anisotropic boost factor within the 60-80 kpc bin in the quadrant polar opposite to the LMC analog at each time step. In addition, these calculations act as proxies for streams in a fully circular orbit centered in each bin in the respective quadrants.
The anisotropic boost factor in m12b both with and without contributed subhalos increases up to 60% as the LMC analog approaches its first pericenter. While both cases for m12b show a positive increase (a positive net inflow of subhalos in the bin), m12b with contributed subhalo trends is consistently higher than m12b without contributed subhalos by approximately 20%. The trend increase in m12b without the contributed subhalos is caused by the analog's effect on the velocity distribution parameters locally. As the contributed subhalos begin to experience tidal stripping, the trends in both cases for m12b converge at around 0.1 Gyr. We note a boost of about 40% right after the first pericentric passage, mostly in Q. ). The mean radial velocity undergoes significant changes in m12b as the LMC analog approaches its first pericenter passage. However, the radial velocity stabilizes again as the LMC analog undergoes tidal stripping. In contrast, m12i exhibits a relatively constant mean radial velocity throughout the time frame (computed within r vir ).
Figure 7. The probability distribution integral from Equation (9) as a function of time computed using the 3D radial velocity distribution (v rad ) in galactocentric coordinates of subhalos in low (10 6 -10 7.5 M e ; red) and high (10 7.5 -10 9 M e ; green) mass bins, and all the subhalos (black) orbiting at a distance of 40-80 kpc from the galactic center. ( ) v rad for different mass bins show similar variation and behavior. An overall increase in ( ) v rad is noted right before the first pericenteric passage at 0 Gyr. The high jitter observed for the 10 7.5 -10 9 M e mass bin arises from small number statistics, which affects the convergence of the integral owing to the low number of subhalos in this high mass bin.
II in the outer bin owing to the reflex motion as the LMC analog moves from Q.III to Q.IV. For m12i (green), the variation in boost factor stays within ±10%, primarily influenced by the halo evolution and the observed anisotropic velocity distribution (E. C. Cunningham et al. 2019).
The LMC analog perturbs the number density n sub of subhalos in two distinct ways.
First, it introduces an anisotropically distributed subhalo population that is dragged along by the satellite's motion. Second, it induces a response in the MW's DM halo. As the LMC approaches its first pericentric passage, the inner galaxy (consisting of the disk and halo system within 30 kpc) reacts more rapidly than the outer galaxy, resulting in a relative displacement between the two reference frames: the inner and outer galaxy frames. This displacement gives rise to a northsouth asymmetry in density, which is referred to as the collective response (N. Figure 9 plots the total number of subhalos within the virial radius of the MW in decades of total subhalo mass for the MW subhalos (black) and MW subhalos combined with LMCcontributed subhalos (green) at equally spaced times: -1 Gyr (left), 0 Gyr (middle), 1 Gyr (right). The insets in each panel plot the fraction of LMC to MW subhalos as a function of each mass decade, and this fraction is found to remain fairly consistent (10%-20%) around 0 Gyr. This indicates that the subhalos contributed by the LMC analog do not have significantly different masses compared to the MW subhalos.
Figure 10 illustrates the variation in the number of subhalos over time in each quadrant for the inner halo (0-50 kpc; top row) and the outer halo (50-100 kpc; bottom row) in the case of m12b with and without contributed subhalos (left and middle columns, respectively), as well as for m12i (right column). The color bar at the bottom represents the quadrant location of the LMC analog in m12b during a time interval of ±0.5 Gyr.
For m12i, the number of subhalos remains relatively consistent across all quadrants over time, except for Quadrant IV in the outer halo, where the systematically increased subhalo population is caused by the presence of an orbiting dwarf galaxy.
In the case of m12b, we identify the specific effects that contribute to perturbations in the number of subhalos within a time interval of approximately ±0.5 Gyr. Q.I has the major contributions from a combination of collective response and leading contributed subhalos.
1. Contributed subhalos. We observe a relative increase in the number of subhalos in the outer halo Q.III (between -0.3 and -0.05 Gyr), Q.IV (between -0.05 and 0.13 Gyr), and Q.I (0.13 Gyr onward) when comparing the left and middle columns. Additionally, there are enhancements in Q.IV of the inner halo at 0 Gyr, corresponding to the satellite being at a pericenter distance of 38 kpc. 2. Collective response. A significant overdensity is observed in the outer halo (middle column), primarily in Q.I and Q. II, as the satellite approaches pericenter in Q.IV (starting from approximately -0.5 Gyr). 3. DM dynamical friction wake. The middle column shows a consistent number of subhalos after the LMC analog leaves a specific quadrant. In the outer halo, the subhalos predominantly trail the LMC analog in Q.III (near 0 Gyr) and Q.IV (between 0.12 and 0.25 Gyr).
These effects collectively contribute to the observed perturbations in the number density of subhalos and their respective quadrants within the specified time interval for m12b.
In this section, we present our results for the stream-subhalo encounter rates for the synthetic streams (Section 4.1) and the MW streams (Section 4.2). These encounter rates are calculated using Equation (8) for each time step along the integrated stream orbits computed by directly integrating Equation (9) without applying any Gaussian assumptions. We subsequently compute the average encounter rates per gigayear using Equation (11). We report the number of encounters per gigayear per stream length ℓ s per maximum impact parameter b max , denoted as N Gyr enc with ˜0 m ¹ . 9 We also utilize the more traditional approach involving setting ˜0 m = , denoted as encounter rates with ˜0 m = , for the anisotropic boost factor as described in Equation (10). For m12b, we further compute the encounter rates by excluding the contributed subhalos, denoted as encounter rates without the LMC. Our simulations 10)) as a function of time for m12b with contributed subhalos (dark blue), m12b without contributed subhalos (magenta), and m12i (green) evaluated in the distance bins and quadrants marked in the top row and in the 60-80 kpc bin in the quadrants polar opposite to the LMC (bottom row). Both cases with and without contributed subhalos exhibit up to a 60% increase in the boost factor around the LMC analog, with peak effect just before T = 0 Gyr. The outermost regions opposite to the LMC analog also experience a higher boost of about 40% after the first passage. m12i, the galaxy with no massive satellites, shows a variation of approximately 10% attributed to halo evolution. These values can be attributed to computations corresponding to a stream assumed to be in a fully circular orbit at the center of each distance bin in the respective quadrant.
are limited by the subhalo resolution and the presence of artificial disruption near the LMC analog (F. C. van den Bosch et al. 2018;S. B. Green et al. 2021), which might lead to a systematic underestimation of the subhalo number densities and therefore the encounter rates. M. Barry et al. (2023) showed that increasing the particle resolution by a factor of 8 did not affect the number density of subhalos significantly.
While impact weighting of the encounter rates-where the contributions of individual encounters are weighted by their expected effects on a stream's morphology-is crucial for predicting observable changes in a stream's structure (D. Erkal & V. Belokurov 2015a; D. Erkal et al. 2016), here we focus solely on the increase in encounters due to the presence of the LMC's subhalos and the MW's response to the LMC. We have shown that ( ) v rad is similar across different mass bins up to a degree of nonconvergence in higher mass bins (Figure 7), indicating that the velocity field is not correlated with the subhalo masses. Additionally, we demonstrated that the LMC analog contributes a consistent fraction of subhalos across all mass decades globally, while the localized contributions can be significantly higher (Figure 9). This highlights that the impact weighting in the presence of the LMC subhalos on a stream morphology would be largely similar to the weighting in the absence of the LMC subhalos.
Moreover, this consistency indicates that our procedure of combining encounter rates across all mass decades is robust. In addition, it mitigates high jitter in the ( ˜)
v rad integral, and given the uniform fractional contribution of subhalos from the LMC analog, it is straightforward to scale these results with n sub for any given mass decade (M. Barry et al. 2023). This approach ensures that our analysis remains valid and reliable, regardless of the specific mass distribution of the subhalos.
We present statistically averaged encounter rates using the 5000 synthetic streams introduced in Section 2.3 as a function of three main properties: (1) galactocentric distance; (2) location on the sky, defined by the quadrant; and (3) angular momentum direction relative to the LMC analog's angular momentum in m12b.
The encounter rates depend linearly on the number density of subhalos n sub and on the velocity dispersion of subhalos in the stream-centric coordinates, both of which decrease as a function of the galactocentric distance. M. Barry et al. (2023) showed that the number density of subhalos is constant as a function of distance from the center out to 50 kpc at the present day; however, we found a decreasing trend in the number density around our integration time ≈5 Gyr before the present day. The top panel in Figure 11 shows the encounter rate per gigayear per stream length ℓ s and impact parameter b max smoothed using the lowess regression (W. S. Cleveland & S. J. Devlin 1988) as a function of a stream's pericentric distance in each quadrant (color-coded) defined in Figure 1, along with a spherically averaged encounter rate without the anisotropic boost ( ˜0 m = ) and the contributed subhalos (gray curve). The bottom panel shows the fractional changes in the rates relative to the spherically averaged rates. The overall number of encounters per gigayear (N enc /Gyr) decreases as a function of distance in all quadrants. The N enc /Gyr values with ˜0 m ¹ (solid lines) consistently exhibit enhancements of 25%-70%, while the rates with ˜0 m = (dashed lines) have a maximum enhancement of only about 25% from the changes in the number density of subhalos.
Q.I and Q.IV curves with ˜0 m ¹ show the highest increment in the rates, with increases of up to 50%-70% (bottom panel). The boost in Q.IV can be attributed to the LMC analog's contributed subhalo increasing the number density of subhalos and changes in the effective g (given by ˜2 m s). Q.I exhibits higher encounter rates due to the collective response and reflex motion. Q.III shows a relatively similar increase at larger distances owing to the contributed subhalos and the influence of the DM wake. Q.II has the lowest encounter rates, as none of the boosting effects are prominent there. On the other hand, m12i (not shown here) has no quadrant dependence, and the encounter rates with m are higher by 20% when compared with ˜0 m = . We note that, even without an LMC analog, ˜0 m = is generally nonzero; similarly, E. C. Cunningham et al. (2019) showed that velocity isotropy changes as a function of position on the sky.
In Figure 12, we plot the fractional enhancements in N enc /Gyr resulting from the presence of the LMC analog's contributed subhalos in m12b. These enhancements are computed by fractional change in the encounter rates with ˜0 m ¹ in the presence of the contributed subhalos and without these subhalos, using the lowess regression. Notably, the trends in fractional enhancements for encounter rates with ˜0 m = show similar patterns. The contributed subhalos have a more pronounced impact on the encounter rates in Q.III and Q.IV, yielding a substantial increase of nearly 15%-40% (Appendix C). This observation aligns with the findings from the first row of Figure 4. In contrast, Q.I and Q.II in the outer halo exhibit relatively lower enhancements, with values not exceeding 12%. Notably, Q.I in the outer halo showcases the least enhancement, and therefore the higher encounter rates seen at larger distances in Q.I in Figure 11 are primarily influenced by the reflex motion effects. We also found a higher localized boost of about 40% for another LMC-analog system in Appendix C, where the analog falls in closer to the present day (M. Barry et al. 2023).
In summary, our analysis highlights the crucial role of subhalo-stream kinematics, particularly the velocity distribution of subhalos, as a key factor in accurately predicting encounter rates. This results in a minimum overall enhancement of 25%. Moreover, we observe a notable dependence in the encounter rates on the location in the sky. Q.IV experiences an additional enhancement of nearly 10%-40%, driven by the presence of contributed subhalos, while Q.I shows increased rates owing to the interplay of the collective response and reflex motion.
We categorize the synthetic streams as prograde or retrograde based on whether their orbital motion is in the same or opposite direction to that of the LMC analog, respectively. Among the 5000 synthetic streams, approximately 50% exhibit retrograde orbits relative to the LMC analog's orbit, while the remaining 50% are in prograde orbits. Figure 13 plots the averaged encounter rates for streams in prograde (green) and retrograde (red) orbits as a function of a stream's pericentric distance in each quadrant (marked in panels in rows 1 and 3). These rates are computed using robust regression (S. S. Wilks 1938; P. J. Huber 1973) in m12b with (solid lines) and without (dashed lines) contributed subhalos. The subpanels below each quadrant (rows 2 and 4) show Δ all as a function of a stream's pericentric distance, which is the difference between the average encounter rates in each quadrant for each orbit type and the spherically averaged encounter rate for both orbit types, with and without contributed subhalos. The magenta bar in specific panels indicates the distance limits of the LMC analog within each quadrant during the calculation of encounter rates. 2) marked with different colors as a function of time within a time interval of ±0.5 Gyr around the first pericenter of the LMC analog. The top row shows the subhalos in the inner halo (0-50 kpc), while the bottom row shows the subhalos in the outer halo (50-100 kpc). In m12b, we see large variations in the number of subhalos as a function of time. For example, prior to pericenter (-0.5 Gyr onward), when the LMC analog is moving from Q.III to Q.IV, the number of subhalos sharply peaks in Q.II (opposite quadrant) in the middle panel. Looking at the left panel, in Q.IV the number of subhalos increases after pericenter, where the LMC analog moves. While in Q.III the number of subhalos decreases after pericenter as the LMC analog leaves Q.III. These changes are the result of the combination of the DM halo response (DM wake and collective response). The lower bar in the m12b panels marks the analog's quadrant information with color in time.
In Q.I, the retrograde orbits with respect to the LMC analog consistently exhibit higher encounter rates than prograde orbits (approximately 0.02-0.08 more encounters or 50% higher than the spherically averaged gray curve shown in Figure 11), regardless of the presence of the contributed subhalos. This is due to the relative bulk motion of subhalos, which aligns in the same direction as the retrograde stream orbits in Q.I, thereby increasing the likelihood of encounters.
In Q.IV, the prograde orbits in the outer regions have a higher number of encounters (approximately by 0.02 Gyr -1 ) owing to their motion aligning with that of the analog, which increases the likelihood of encountering subhalos, as the streams spend more time in close proximity to the contributed subhalos. Additionally, the presence of a transient DM wake can further contribute to the increased encounter rates for prograde orbits in the outer regions of Q.III.
In Q.II, retrograde streams have systematically lower rates than both the spherically averaged rates and prograde orbit rates. We think that this may be due to the fact that streams move between quadrants over the full integration time of 0.5 Gyr. The retrograde streams in Q.II preferentially enter via Q.III, where the relative velocity with subhalos from the LMC analog is high since most subhalos in that quadrant have a bulk velocity in the same sense as the analog's orbit. This then produces an apparent decrement in retrograde encounters in Q. II. Conversely, streams that are moving prograde in Q.II at the end of the integration period came from Q.I, which is not strongly affected by the analog's contributed subhalos. This would also produce the difference between the rates calculated with and without including the contributed subhalos that we also observe in Q.II. This quadrant serves as a reminder that the impact of the LMC on stream-subhalo encounters is not quite as simple as we would hope! m ¹ (solid) and ˜0 m = (dashed). The rates are plotted as a function of a stream's pericenter distance at T = 0 Gyr in each quadrant (marked by different colors). The gray curve marks the encounter rates for m12b ˜0 m = and without the contributed subhalos, as well as for m12i (rescaled to match the m12b halo mass), averaged over all quadrants. These curves are estimated using a lowess regression (W. S. Cleveland & S. J. Devlin 1988) from the 5000 synthetic streams. Streams in Q.I and Q.IV have the highest probability of encounters, as the LMC analog's position in m12b at T = 0 Gyr leads to a higher concentration of subhalos in Q.IV (at large distances), and the majority of the reflex motion and collective response effects are dominant in Q.I. Outer regions of Q.II have the lowest encounter rates of all, consistent with the spherically averaged rates marked by the gray curve. The second row plots the fractional changes in each encounter rate curve with respect to the spherically averaged rates.
Figure 12. Fractional enhancements in the encounter rates due to the contributed subhalos from the LMC analog as a function of a stream's pericentric distance in each quadrant (marked with different colors) for m12b. Streams exhibit a consistent increase of approximately 20% in encounter rates across all distance ranges in Q.III and inner regions of Q.IV owing to the contributed subhalos. In contrast, outer regions of Q.I show a smaller increment of ∼10% owing to contributed subhalos. Notably, Q.II has a constant increment of about 12% at larger distances, indicative of the presence of trailing contributed subhalos. These subhalos can boost the encounter rates through the halo by 15%.
Figure 13. The average encounter rate (rows 1 and 3) is estimated using separate robust regressions (S. S. Wilks 1938; P. J. Huber 1973) for prograde (green) and retrograde (red) streams, taking into account the stream's pericenter distance and sky position (quadrants marked on the panels). These calculations are performed in m12b with (solid) and without (dashed) contributed subhalos from the satellite. The difference with respect to the spherical average, denoted as Δ all , is also shown in rows 2 and 4 in each quadrant for each curve. The magenta bar marks the LMC analog's distance in each quadrant within the time frame of the encounter rate calculation. Additionally, predicted encounter rates for the MW streams GD-1, Orphan, and AliqaUma, categorized by their orbit types, are also shown in their respective quadrants, plotted as a function of their present-day distance. We note that streams in retrograde orbits with respect to the LMC analog consistently exhibit higher encounter rates in Q.I (regions polar opposite to the analog's trajectory).
Furthermore, the encounter rates with ˜0 m = (not depicted) show no discernible trends between different orbit types across any of the quadrants. However, the introduction of a nonzero m complicates the systematic perspective owing to the influence of the stream's apparent motion on the anisotropic boost factor. Additionally, there is no observed dependence or correlation between a stream's orbital eccentricity and encounter rates.
We annotate the GD-1 stream (C. J. Grillmair & O. Dionatos 2006), the Orphan stream (C. J. Grillmair 2006;V. Belokurov et al. 2007;S. E. Koposov et al. 2019), and the AliqaUma stream (T. S. Li et al. 2021b) in their respective quadrants and orbit types relative to the LMC analog. The trends observed in our synthetic streams generally align with those observed in the MW streams. The Orphan-Chenab stream closely interacts with the LMC and extends across both hemispheres of the MW (S. E. Koposov et al. 2023). Additionally, S. Lilleengen et al. (2023) have demonstrated the effects of the deforming MW halo on the proper motion and morphological structure of streams, making it challenging to differentiate between changes caused by subhalo-stream interactions and other influences.
We have observed that synthetic streams around halos with an LMC analog display a diverse range of encounter rates, contingent on several factors: (1) the specific value of m associated with the stream; (2) the stream location in the sky; and (3) the stream orbital rotation direction in relation to the LMC analog's orbit. Notably, our findings demonstrated boosted encounter rates for orbits located in the northern hemisphere and in regions near and diametrically opposite to the LMC.
Expanding on our analysis, we now proceed to calculate encounter rates per gigayear per ℓ s and b max for all the MW streams within our simulation. We employ the same orbital integration methodology as employed for synthetic streams. However, for the MW streams, we substitute the initial conditions with the present-day phase-space coordinates of each stream in the MW. To establish these initial conditions, we utilize the median position and velocity of each stream track. We further adjust these conditions by applying a rotation that emulates the orientation of the actual streams relative to the real LMC as discussed in Section 2.3.
Our results are summarized in Table 3; we present the computed encounter rates (N enc /Gyr) in both m12b and m12i for the known MW streams (C. Mateu 2023). The table also includes the orbital properties of each stream, such as its present-day distance from the host galaxy (d) and the sky quadrant (Q) it occupies at T = 0 Gyr in the context of m12b. The eccentricity of the orbit (e) and its grade (G) with respect to the LMC analog in m12b are also provided. The grade indicates whether the stream's motion is prograde (in the same direction as the LMC analog's orbit, denoted by "+") or retrograde (opposite to the LMC analog's orbit, denoted by "-"). These properties are calculated using the m12b potential models of DM halo and baryonic disks from A. Arora et al. (2022). It is crucial to emphasize that around 75% of the MW streams are situated within a radius of 20 kpc at T = 0 Gyr. While interpreting these encounter rates, it is important to be cautious owing to the limited resolution in the inner regions. M. Barry et al. (2023) showed that the total subhalo population within a 50 kpc radius, 5 Gyr ago, was approximately twice as high for subhalos with M sub 10 6 M e compared to the present-day population. As a result, we can adjust our encounter rates by dividing them by a factor of 2. For example, M. Barry et al. (2023) estimated that the GD-1 stream, with a length of 15 kpc and an average impact parameter of 1.5 kpc, could experience approximately five encounters per gigayear using traditional symmetry and isotropy assumptions. With the inclusion of the LMC-induced subhalos, this rate increases twofold. In our analysis, our rescaled factor for present-day encounter rate for GD-1 with the same stream length and impact parameter is 3.5 encounters per gigayear for m12i and 4.3 encounters per gigayear for m12b. These numbers are more or less similar because m12i has about 1.2 times fewer subhalos than m12b without the LMC-contributed subhalos. Given the location of GD-1 in the sky and distance from the center, we don't expect any boosts in encounter rates due to the LMC, and our numbers are consistent with M. Barry et al. (2023) to the first order.
Among the MW streams analyzed, the Aquarius-W11 stream (M. E. K. Williams et al. 2011), Lethe-G09 (C. J. Grillmair 2009), and Phlegethon-I21 (R. Ibata et al. 2021) in the inner halo (within 30 kpc) exhibit the highest encounter rates within the m12b simulation. Conversely, the Cetus-Y13 stream (W. Yam et al. 2013) has the highest encounter rate in the outer halo (at distances greater than 30 kpc). The streams C-7-I21 (R. Ibata et al. 2021), Ophiuchus-C20 (N. Caldwell et al. 2020), and OmegaCen-I21 (R. Ibata et al. 2021) have the lowest encounter rates among the analyzed streams. These streams are located in the inner halo, orbiting within 7 kpc of the galactic center. Additionally, the Eridanus-M17 stream (G. C. Myeong et al. 2017), positioned at a distance of 100 kpc, exhibits the lowest encounter rates compared to other streams.
Overall, we observe consistent encounter rate trends for the MW streams across m12i and m12b, both with and without m. Streams situated beyond a present-day distance of 20 kpc in the m12b exhibit a clear dependence on location, with enhanced encounter rates in Q.II and Q.IV. Conversely, in m12i streams demonstrate no azimuthal dependence in the encounter rates. Incorporating m yields approximately a 30% boost for both m12b and m12i. The sense of rotation of the streams does not notably impact the encounter rates, mainly due to the proximity of these streams to the galactic center, while these effects are prominent for streams positioned in the outer halo regions.
Currently, there are ≈100 observed streams in the MW. However, in the next decade we expect to discover more streams, especially in the outer halo, thanks to upcoming surveys such as the LSST (Ž. Ivezić et al. 2019). Furthermore, we will have a multidimensional view of the phase space and kinematics of the streams with unprecedented accuracy. As such, the search for disturbances from DM subhalos can be done systematically across large areas of the sky. With the insight from previous sections, we can anticipate specific areas in the MW sky where the likelihood of finding streams perturbed by subhalos is higher.
In Figure 14, we highlight the specific regions in the sky with an increased probability of detecting stream-subhalo interactions. These regions are identified by analyzing the sky positions of synthetic streams at T = 0 Gyr in m12b, rotating them such that the LMC analog matches the real LMC location in the MW, considering the expected changes in the rates with and without the contributed subhalos. These regions are shaded manually by running a k-means clustering algorithm on the stream locations (J. MacQueen 1967), dividing the stream positions in the Galactocentric coordinates into three distinct groups: one group for the background, another for the reflex motion, and a third for the LMC subhalos and DM wake categories, while also removing any outliers.
The red shaded region is the polar opposite region to the LMC, where the dominant effect on stream-subhalo interactions is the reflex motion response, in both position and velocity. The purple shaded region indicates a higher subhalo density due to the presence of LMC subhalos in the leading debris and the LMC itself. Streams located within these shaded regions are expected to have a higher likelihood of encounters with subhalos. Our estimates based on Section 4 suggest an enhancement ranging from 50% to 70% (depending on the stream's orbital orientation; Figure 13) in the reflex motion region (red shaded area) and from 20% to 50% in the region influenced by the leading LMC subhalos (see Figures 11 and 12) (purple shaded area) with respect to the spherically averaged rates without accounting for the LMC-contributed subhalos.
Additionally, we include the MW streams orbiting beyond 20 kpc from the galactic center, with encounter rates greater than the 90th percentile limit (4.4 encounters per gigayear) based on our calculations in m12b (see Table 3). We also list their corrected-for-present-day (reduced by a factor of 2) encounter rates per stream length of 10 kpc and impact parameter.
We predict that streams such as NGC1261-I21 (R. It is important to note that caution is required when interpreting these values as calibration benchmarks for specific stream regions. Instead, they offer insight into the potential systematic uncertainties that can arise from departures from equilibrium assumptions in these measurements.
In this paper, we have investigated the impact of a massive satellite on the encounter rates between stellar streams and DM subhalos. Leveraging two representative systems from the FIRE-2 hydrodynamical zoom-in simulations-one featuring an LMC analog and the other lacking any massive satellite-we address two key questions. First, we evaluated how the LMC analog's subhalo population contributes to enhancing encounter rates with host stellar streams, building on the work of M. Barry et al. (2023) and using the same simulations. Second, we investigated how the host halo's response influences encounter rates galaxy-wide.
Our findings indicate a general decrease in encounter rates with increasing pericentric distance, independent of the eccentricity of a stream's orbit. The presence of a massive satellite, such as the LMC analog, introduces an anisotropic boost, resulting in varying encounter rates for streams across different regions of the sky relative to the satellite. Key factors influencing these asymmetric effects are ranked as follows:
1. Mean radial motion. The mean radial motion of subhalos with respect to streams, represented by m, consistently enhances encounter rates by at least 30% even in systems without the LMC analog. Under the influence of an LMC analog, localized boosts can jump up to 70% as observed in Q.I, compared to the spherically symmetric and ˜0 m = rates (Figure 11). 2. Contribution from the LMC-analog subhalos. The LMC analog brings its own set of subhalos with unique phasespace orbits, as depicted in Figure 4. These subhalos can Gaia et al. 2018). Shaded regions on the map indicate higher overall probabilities of encounter rates, which have been determined using the top 1% synthetic streams in m12b, rotated to align with the MW-LMC orientation at the present day. To aid visualization, the synthetic streams have been color-coded using a k-means clustering algorithm (J. MacQueen 1967), and shaded regions have been manually drawn. Also marked are the MW streams with the highest encounter rate enhancements located at distances greater than 20 kpc, along with their reduced (by a factor of 2) rates per 10 kpc stream length and impact parameter in m12b. Streams within the shaded regions labeled as reflex motion (red) and subhalos brought in by the LMC (purple) are anticipated to have a higher probability of subhalo-stream interactions by as much as 50%-70% and 20%-50%, respectively, compared to the spherically averaged rates without accounting for the LMC-contributed subhalos.
boost the encounter rate up to 10%-40% in Q.III and Q. IV, aligned with the LMC analog's motion (Figures 12 and 17). This boost may be much higher in the MW today (Appendix C). M. Barry et al. (2023) showed that the number density can be boosted up to a factor of 2 owing to the LMC-mass satellites. 3. Orbital alignment. Streams in retrograde orbits with respect to the LMC analog, particularly located in the outer halo and the opposite hemisphere from the analog, exhibit enhanced rates up to 30%, with the most substantial boosts occurring in Q.I compared to prograde streams (Figure 13).
In conclusion, our investigations reveal that streams situated in proximity to the LMC and in regions opposite to it are more likely to interact with subhalos, potentially leading to detectable morphological changes. Building on this insight, we have identified regions within the MW-LMC celestial sphere likely to display signatures of stream-subhalo interactions, as depicted in Figure 14. Detecting such signatures and constraining the perturbing subhalo masses will offer valuable insights into allowed DM models with specific subhalo mass scales. Furthermore, our results emphasize the complexities introduced by the presence of a massive satellite such as the LMC, suggesting that the MW might not be an ideal laboratory for constraining the subhalo mass function, due to the need to consider departures from the equilibrium assumptions. In this context, finding extragalactic streams orbiting in galaxies in a state of dynamical equilibrium presents promising opportunities for precise constraints on the DM subhalo mass function. With upcoming missions like the Roman Space Telescope (D. Spergel et al. 2013Spergel et al. , 2015)), we anticipate the detection of such streams in other galaxies (S. Pearson et al. 2022;C. Aganze et al. 2023), enhancing our ability to derive valuable insights into the nature of DM. moving away from the stream) and diverges to infinity for a large negative mean velocity (all the subhalos approaching the stream in the finite time interval). In the MW, the velocity dispersion of the outer halo is of order σ = 100 km s -1 (A. J. Deason et al. 2012;J. G. Cohen et al. 2017), and the expected mean velocity is ˜35 m = km s -1 (D. Erkal et al. 2019; M. S. Petersen & J. Peñarrubia 2020). For an MW stream in a circular orbit the anisotropic boost factor will scale the encounter rates to 1.75 or 0.5 depending on whether m is locally negative or positive in the stream's frame. Notably, a negative m value leads to a substantial rate enhancement, and a positive m leads to a reduction. For instance, a g of -0.25 results in a 50% boost in encounter rates, while a + 0.25 value leads to a 35% decrease. For halos with no massive satellites, M. Barry et al. (2023) showed that in the galactocentric frame at the present day the orbital velocity vectors of subhalos are generally isotropic at all distances.
We examine the validity of assuming a Gaussian distribution ( ˜˜) J , gauss m s for ( ) P v R . It is noteworthy that the pdfs subjected to analysis in this section exclusively concern subhalos positioned around streams, in accordance with the description provided in Section 2.5 within the context of the m12b at T = 0 Gyr. This specific choice is motivated by the fact that the LMC analog's first pericentric passage induces significant perturbations on the Gaussianity assumption, due to its pronounced dipolar moment (E. C. Cunningham et al. 2020; M. S. Petersen & J. Peñarrubia 2020).
The left panel of Figure 16 shows ( ) P v R (solid lines) and their corresponding Gaussian fits (dashed lines) for a representative set of MW streams (varied colors). Notably, the means of different pdfs and their fits display similar trends. It is important to note that this mean value can shift, toward either more negative or more positive values, thereby leading to the conditions f (γ) 1 or f (γ) 1, respectively. This shift's potential influence on the scaling factor behavior is highlighted. The alignment between Gaussian distributions and actual pdfs visually confirms the Gaussian fitting assumption.
For a quantitative analysis, we perform numerical integration for the 5000 synthetic streams at T = 0 Gyr in m12b, assuming a Gaussian approximation (see Equation (10)), hereafter labeled as gauss , and conducting the exact numerical integral
) without any assumption (see Equation (7) in Section 2.5). We compute the residual as 1 gauss num -to assess the errors introduced by the approximation.
The bottom panel of Figure 16 shows the resulting residual distribution for the 5000 synthetic streams in m12b at T = 0 Gyr. The dashed line represents the 50th quantile, and the dashed-dotted lines mark the 16th and 84th quantiles. The majority of errors fall within 10% of the expected values, with the mean and standard deviation of the distribution being 0.01 and 0.08, respectively. This translates to an average error of 1% with a deviation of 8%. We also note that the majority of higher errors (10%) occur for streams close to the galactic center and the disk stemming from the low subhalo population. These results highlight the accuracy of the Gaussian approximation method, especially compared to other assumptions used in such analysis.
For a qualitative evaluation, we test the Gaussian assumption using the Shapiro-Wilk test (S. Remarkably, the p-value counts demonstrate a logarithmic decrease, highlighting a gradual transition from non-Gaussian to Gaussian behavior in the pdfs. This observation further endorses the Gaussian approximation's appropriateness for the pdfs. While the numerical 1D integrals are efficient, the assumption of Gaussianity remains valid and imposes no significant constraints compared to other approximations commonly employed in such analyses. f g in Equation (10) that is dependent on the Gaussian parameters of the cylindrical radial velocity distribution of subhalos with respect to a stream as a function of a typical g , i.e., is a function of mean radial velocity ( m) and velocity dispersion (s) of subhalos in the stream-centric coordinates. Theoretically the factor converges to 0 with ˜0 g and diverges to infinity for ˜0 g . Parameter ( ˜) f g strongly depends on m and s, and hence assuming ˜0 m = can bias the estimates of encounter rates.
http://fire.northwestern.edu/latte
The Astrophysical Journal, 974:286 (22pp), 2024 October 20 Arora et al.
The notation ˜0 m ¹ is used here, but the rates are computed by directly integrating Equation (9) without applying any Gaussian assumption.