The Resolution Paradox: Navigating the Inverse Problem of Dark Matter Subhalo Detection

Abstract

We develop a custom N-body simulation to quantify the relationship between dark matter subhaloes and gaps in stellar streams. Using a Java-based gravitational simulator with semi-implicit Euler integration and cubic softening, we model interactions across parameter spaces of 5-50 stream particles and 0-10 subhaloes (0-5 for the 5 particle case). After detecting gaps and predicting subhaloes through Python analysis, statistical analyses in R reveal a strong positive correlation between subhalo abundance and gap counts (), though systematic overprediction of subhalo density occurs at low subhalo abundance. Counterintuitively, prediction error in subhalo counts increases with stream resolution despite stronger correlations. These findings demonstrate that stream gaps reliably indicate subhalo presence but require more sophisticated inference methods beyond mere gap counting; nevertheless, this study provides a foundation for detecting dark matter substructure via stellar streams and advancing our understanding of hierarchical galaxy formation.

Keywords: Dark matter subhaloes; stellar streams; tidal disruption; gravitational perturbation

Introduction

Dark Matter (DM) dominates the matter content of the universe, constituting approximately 84.4% of all matter¹; its gravitational influence results in the formation of cosmic structures ranging from galaxies to the smallest subhaloes. Cold Dark Matter (CDM) simulations predict a hierarchy of DM substructures within galactic halos^2‘3. These subhaloes, though invisible, leave imprints on visible tracers, notably stellar streams, and on larger scales, gravitational lensing. While lensing of background galaxies reveals DM distributions via spacetime curvature⁴, it struggles to resolve subhaloes below⁵. Stellar streams, by contrast, respond to these small scales directly.

Stellar streams, essentially elongated filaments of stars formed from the tidal disruption of globular clusters as well as dwarf galaxies, act as sensitive gravitational probes (see Figure 1). Because they have a low velocity dispersion, they are ideal for detecting more subtle perturbations caused by DM subhaloes⁶.

Stellar streams interact with dark matter substructure through two distinct but interconnected gravitational processes. First, direct gravitational perturbations occur when subhaloes pass near stream stars and deliver impulsive velocity kicks that displace stars from their original orbits⁸. Second, tidal forces arise from the differential gravity across the stream’s width and stretch the stellar distribution into varying morphologies⁷. Perturbations dominate during close encounters, but tidal effects become increasingly important for the stream’s long-term evolution, particularly when considering: (1) the stream’s extended geometry, and (2) repeated interactions over time^9‘10. This duality is apparent in the Pal 5 stream, which shows sharp density gaps from recent encounters and broad fanning from cumulative tidal distortion. Gaps and spurs in streams such as GD-1 and Pal 5 are thought to arise from gravitational interactions with subhaloes⁸, but although cosmic simulations (e.g., Aquarius) model these types of interactions, they often require immense computational resources and also conflate baryonic effects with DM dynamics¹¹.

Still, even with advances in these simulations, few tools exist to invert observed stream perturbations into constraints on subhalo properties. Current methods often rely on simulating subhalo impacts and comparing them to observations but lack a systematic approach to solve the inverse problem¹². To address the limitations of current approaches, a minimalist Newtonian simulator that isolates the gravitational interaction between subhaloes and stellar streams is developed in this study—appropriate given the non-relativistic velocities and spatial scales that are typical of Milky Way stellar streams¹³. Indeed, General Relativity plays a significant role in cosmological dynamics; however, Newtonian gravity suffices for modeling local perturbations.

This simulation is built on two core ideas from dark matter theory.

Hierarchical structure formation suggests that subhaloes are leftover fragments from a long history of mergers, as predicted by cold dark matter (CDM) models. The initial positions assigned to them in the code (sx,sy,sz) are meant to reflect this accretion history¹⁴.

Tidal disruption theory explains how stellar streams form: when tidal forces from the host halo become strong enough to overcome a satellite’s self-gravity, stars are stripped away and form elongated streams¹⁵.

This model takes into consideration the effects of tidal disruption and gravitational dynamics, all while accounting for forces from both the host galaxy and passing subhaloes of varying proximities and velocities. By omitting baryonic physics and relativistic effects, this model isolates subhalo properties’ influences on stream morphology and, in turn, provides a more focused view on the inverse problem. This connection between theoretical subhalo populations and observational stream features could potentially refine estimates of the Milky Way’s DM substructure abundance and directly test CDM predictions¹⁶.

This study pursues three primary goals:

1. Forward Modeling: Simulate the effects of subhalo flybys on stellar streams under varying parameters (e.g., subhaloes and particle number) in the Milky Way.
2. Inverse Problem Exploration: Investigate whether gap size within streams can be used to infer the density of perturbing subhaloes.
3. Model Validation: Qualitatively evaluate whether this simulation outputs align with physical expectations.

Accordingly, guiding this work is the central research question: How can we model the perturbative effects of dark matter subhaloes on stellar stream orbits using Newtonian dynamics in order to infer their distribution in galaxies?

Besides its scientific contributions, this project also has pedagogical value; the modular Java classes used in this study demonstrate central N-body techniques in an accessible way that is suitable for classroom adaptation¹⁷.

This simulation is also intentionally simplified—modeling with symmetrical and collisionless subhaloes, a static host potential, and an omission of baryonic physics and relativistic effects—to isolate the essential gravitational interactions between subhaloes and stellar streams. Furthermore, the particle and subhalo counts are lower than those in the actual Milky Way to streamline computation, and the symmetry for each particle cluster and subhalo reduces complexity without obscuring this experiment’s main objectives. These simplifications are all enforced deliberately to create a balance between run-time efficiency as well as this model’s physical interpretability, though future extensions could incorporate more accurate numerical schemes such as leapfrog integration or tree-based force solvers to improve model precision¹⁸.

Working within these simplifications, the research involves a numerical simulation, perturbation analysis, and statistical validation. A custom N-body simulator is first developed in Java to model gravitational interactions between subhaloes and stream particles across varied parameter spaces. The resulting simulation data is then processed through a Python-based analysis that utilizes density-based gap detection and velocity kick identification algorithms. Finally, comprehensive statistical analyses are conducted in  to quantify correlations between actual subhalo abundance and detected stream perturbations and to evaluate the efficacy of the inverse problem approach for subhalo detection.

Methodology

Introduction to Computational Approach

This study used a custom Newtonian N-body simulation framework consisting of two integrated components: a Java-based gravitational dynamics simulator that modeled the interactions between dark matter subhaloes and stellar streams, and a Python-based analysis pipeline that detected gravitational perturbations and predicted subhalo abundances from stream morphology. The Java simulation implemented a softened Newtonian gravitational potential, where  simulation units represented the gravitational softening length, which was chosen to prevent extreme numerical divergences and to produce more interpretable graphs.

Gravitational Softening

The gravitational force between two particles (or subhaloes) was computed using a softened inverse-square law^19‘20:

   

Here, G was the gravitational constant (set to 1 for simplicity), and were the masses of the interacting objects in simulation units, and and were their position vectors. was the softening length. After a series of calibration tests for multiple particles and subhaloes, it was found that smaller values (eg,. ) resulted in orbits extremely prone to numerical instability during close encounters, while excessively large values (e.g., = 3) overly suppressed gravitational interactions. The chosen value of 1.75 was the smallest softening value to guarantee numerical stability over the full simulation runtime. This softening law was applied to a three-dimensional space instead of a two-dimensional space; a similar transition from stable elliptical orbits to rosette-style patterns as softening increased can be seen in Figure 2, which presents both cubic and Plummer softening formulations.

Indeed, Plummer softening is more standard, but the cubic form was selected for its computational efficiency in avoiding a square root operation²¹. This resulted in a stronger suppression of the force at small separations () compared to Plummer softening, which further stabilized close encounters. But overall, the large-scale dynamics involved in gap formation remained unchanged. Figure 2 demonstrates a minor change from one method to the other: cubic softening produced clean elliptical orbits characteristic of Keplerian motion, whereas Plummer softening introduced slight orbital precession due to its different small-scale force behavior.

Time Integration

To evolve the system forward in time, this simulation used a semi-implicit Euler method (also known as the symplectic Euler method), which offered better energy conservation than the standard explicit Euler scheme for gravitational systems. For instance, the Forward Euler method accumulated exponential error over time, leading to instability in long-term simulations. In contrast, the semi-implicit Euler method helped mitigate these issues by updating the velocity using the current acceleration, which provided greater numerical stability. The Backward Euler method’s system of equations, while more stable, was computationally expensive when solved iteratively through methods like Newton’s method. The update rules in the semi-implicit method for position position r and velocity v were:

   

   

where was the acceleration at timestep n, and was the fixed time step. The second equation updated the position r for the next time step; appeared on the right-hand side, which indicated a form of implicit updating for the velocity²².

Though higher-order methods such as leapfrog integration are standard for cosmological simulations, their superior long-term energy conservation was not a primary requirement for this study’s short-duration encounters. To validate the semi-implicit Euler approach, a leapfrog integrator, implemented for comparison in Figure 3, confirmed that both methods would produce almost identical results for key outputs: gap counts (11 for semi-implicit, 12 for leap-frog) and predicted subhaloes (4 for both); additionally, there was only a slight variation in terms of velocity kicks. Given that leapfrog integration yielded similar results, other energy-conserving methods such as Runge-Kutta would also be expected to perform as such²³. Therefore, the semi-implicit Euler method was kept for its implementational and conceptual simplicity, ultimately aligning with the pedagogical goals of this study without sacrificing the accuracy required for analysis.

Tidal Perturbation

The tidal acceleration exerted by a subhalo on a stream particle was derived from the tidal tensor formalism²⁴. For a subhalo of mass at a distance r, the tidal acceleration in the direction of the stream’s angular momentum was:

   

This term was computed at each timestep and applied as a perturbation to the particle’s velocity, essentially modeling the differential gravitational stretching that leads to gap formation.

Stellar Stream Initialization

Stellar streams were modeled as collections of test particles (5-50 particles) initialized in relatively close proximity to enable collective stream-like behavior; each particle represented a cluster of stars. The upper limit of 50 particles was chosen, as it was sufficient to capture the bulk dynamics and tidal stretching while avoiding further numerical instabilities that appeared with including hundreds of particles-which also would have complicated our goal of isolating the fundamental relationship between subhaloes and stream perturbations. Hence, this approach prioritized a clear understanding of the mechanism over simulating a fully realistic stream.

Each particle mass was set to approximately 300 times less than subhalo masses (0.1 vs. 30 simulation units) to make sure the stream would respond to subhalo perturbations. Initial velocities were deliberately set with higher magnitudes than subhalo velocities ( simulation units compared to subhalo ) such that the stream particles would maintain orbital motion around the galactic center while remaining susceptible to gravitational perturbations from passing subhaloes. This velocity configuration, with minimal x and z components (), resulted in primarily planar initial orbits with slight vertical motion to enable three-dimensional interaction effects.

DM Subhalo Initialization

Subhalo masses were set to approximately 6% of the galactic mass (30 vs. 500 simulation units) to produce clearly observable perturbations. For a Milky Way-like host halo mass of , this ratio corresponded to subhalo masses of , which placed these subhaloes at the high-mass end of the substructure spectrum. The resulting cumulative mass fraction in substructure was therefore 6%-higher than the typical 1% predicted for Milky Way-mass haloes but fell within the range of 5-10% established for group-to-cluster mass haloes in cosmological simulations²⁵. Although this subhalo and galaxy mass ratio was higher than that of typical ΛCDM subhaloes, it produced a stronger signal-to-noise ratio for gap formation, which provided a foundational understanding of the gap formation mechanism that can be scaled to the more cosmologically relevant mass range of of  in future work²⁶. Subhaloes were positioned initially outside the stellar stream orbits to maximize gap creation through gravitational encounters as opposed to capture or binding effects. Velocity vectors were randomized in direction with alternating orbital inclinations (positive and negative y and z components) to simulate a diverse population of subhalo orbits; simultaneously, lower velocity magnitudes than stream particles were maintained () to ensure sufficient interaction timescales for gap formation. This idealized velocity set-up aimed to further isolate the physics of tidal and gravitational forces, and was independent of the often unknown and extremely variable initial conditions of real streams.

Data Collection on Java

The Java simulation advanced through discrete time steps, where each step represented a timestep of in this simulation’s time units. The duration of 5.0 simulation time units allowed the system to capture multiple close encounters between subhaloes and the stellar stream. The Java simulation state was then recorded to a CSV file at regular intervals of every time steps; the file captured the precise position (x,y,z), velocity (vx,vy,vz), angular momentum components (lx,ly,lz), their derivatives (lvx,lvy,lvz), and mass for every particle and subhalo at each snapshot. The number of test particles in the stream was varied across the values 5, 10, 15, 20, 30, 40, 50. For each particle count, the number of actual subhaloes was varied from 0 to 10. However, for the 5 particle resolution, the subhalo range was limited to 0-5 due to numerical stability concerns and to maintain physical plausibility-specifically, to make sure subhaloes did not outnumber stream tracer particles, which would violate the assumption that streams serve as sensitive detectors of substructure.

Each CSV file was run through Python to contribute to another data file, categorizing particles, gap_counts, velocity_kicks, predicted_subhaloes, actual_subhaloes, predictedsubhalodensity_actualsubhalodensity_ratio; every combination of subhalo and particle values was tested three times, and the average value for each category per combination was inputted into this data file. An additional two columns would be added on the data sheet, one for standard deviation per three trials for every combination of subhalo and particle values, as well as the standard mean error. Lastly, the data file was read by R to run statistical analyses.

Gap and Subhalo Detection Algorithm

The Python analysis utilized a density-based gap detection algorithm that first computed stream density profiles using linear kernel density estimation along the stream’s longitudinal axis with 50 bins. Regions where density fell below 50% of the mean stream density were flagged as gaps, representing potential subhalo-induced perturbations²⁷. This 50% threshold was determined after sensitivity trials revealed it to be an optimal compromise. Thresholds lower than 30% were too strict and missed a large amount of perturbations, leading to a severe underestimation of subhaloes; meanwhile, higher thresholds above 60% were too sensitive to random fluctuations and caused a severe overestimation.

Velocity kicks were detected from velocity changes ( arbitrary simulation units) between timesteps solely to provide additional evidence of gravitational encounters. The cutoff was also determined through a sensitivity analysis. A threshold that was too low () failed to filter out the stream’s intrinsic velocity dispersion, which yielded a high background of detected kicks even in unperturbed streams. Conversely, a threshold that was too high () was overly restrictive and discarded subhalo-induced kicks, thus providing insufficient evidence of gravitational interactions. The chosen value of 2.0 resided within an optimal range where the kick count meaningfully correlated with subhalo encounters.

Subhalo abundance was predicted using a calibrated linear relationship: Predicted Subhaloes = Gap Count / k, where the calibration factor k was determined empirically. To establish this factor, we analyzed the high-signal regime of our parameter space-simulations with the maximum of 10 subhaloes (5 subhaloes for the 5 particle resolution). Across all stream resolutions, these simulations consistently produced a total gap count that was, on average, three times the number of actual subhaloes. This resulted in a calibration factor of k = 3.

Density and Statistical Analysis

The spherical shell volume for density calculations was defined with a radius 20% larger than the maximum subhalo galactocentric distance, with a thickness equal to 5% of this radius. Density ratios were computed as the ratio of predicted-to-actual subhalo densities within this volume-unique for each Java simulation trial. Statistical analyses in R included Pearson correlation tests between actual subhalo counts and gap detections across stream resolutions, t-tests distinguishing gap counts in perturbed versus unperturbed streams (p 0.05 threshold), and linear regression quantifying the gap-subhalo relationship. Prediction accuracy was assessed through mean error analysis and density ratio comparisons, with all analyses accounting for the experimental design of 3 runs per parameter combination across 7 different particle resolutions.

Results

To establish a baseline for the gap detection method, the algorithm was first applied to a control simulation of an unperturbed stellar stream. As shown in Figure 4, a stream was initialized with  test particles in a cold configuration, orbiting the galaxy center in the complete absence of dark matter subhaloes. Any gaps detected in this scenario signified false positives. This false positive rate must be accounted for when solving the inverse problem of inferring subhalo abundance from gap counts in perturbed streams. Velocity kicks were determined from numerical integration and time sampling, and the resulting distances were converted to kiloparsecs (kpc) using a scaling factor of 3. This established a realistic orbit with a galactocentric radius ranging from approximately10 to 40 kpc²⁴.

To test the algorithm’s sensitivity, a single perturbing subhalo with a mass equivalent to 6% of the galaxy center was introduced to the 5 particle stream configuration. The introduction of a new element produced visual and measurable perturbations: 2 additional gaps and 4 additional velocity kicks compared to the unperturbed baseline case. Despite the signs of an interaction, identical integer subhalo estimates were predicted for both the perturbed and unperturbed case (see Figure 5).

Despite this limitation at low subhalo abundances, a relationship between subhalo abundance and gap formation was still able to be quantified across a range of stream resolutions. As shown in Figure 6, actual subhalo counts and detected gaps showed a significant positive correlation (), suggesting that gap formation scales systematically with subhalo abundance.

The strength of the correlation between subhalo abundance and detected gaps was analyzed as a function of stream particle resolution. All tested resolutions yielded correlations with p-values below 0.05, and r values ranging from 0.759 to 0.914. The variation in correlation strength revealed an optimal detection range around 40 particles (); lower particle counts (5-20) showed slightly more variable correlations, while higher resolutions (30-50 particles) maintained relatively strong detection performances (see Table 1).

Particles	Correlation (r)	P-value	Sample Size (N)
5	0.893	0.016	18
10	0.895	0.00020	33
15	0.831	0.0015	33
20	0.759	0.0067	33
30	0.860	0.00068	33
40	0.914	0.000083	33
50	0.870	0.00049	33

However, analysis of the predicted subhalo densities implied a systematic bias-namely, a consistent overprediction (ratio 1) for lower subhalo counts (1-7), with particularly severe overestimation for single subhaloes (ratio 5.6). There was only a progressive convergence toward accurate prediction (ratio 1) as well as a progressive decrease in standard error at high subhalo abundances (9-10), as detailed in Table 2. All ratios were less than p = 0.05 different from perfect prediction (ratio 1) except for the 10 subhaloes case where p 0.082.

Actual Subhaloes	Mean Density Ratio	Standard Error	p-value	Sample Size (N)
1	5.6	0.528	0.00013	21
2	3.0	0.289	0.00045	21
3	2.0	0.138	0.00040	21
4	1.7	0.130	0.0020	21
5	1.4	0.081	0.0018	21
6	1.9	0.083	0.000091	18
7	1.8	0.088	0.00025	18
8	1.6	0.070	0.00046	18
9	1.2	0.089	0.034	18
10	1.1	0.061	0.082	18

A further unexpected finding concerned stream resolution. High resolution (30-50 particles) tended to improve correlation strength in terms of statistical significance, but also led to an increase in mean absolute prediction error, the highest mean error being 4.72 at 40 particles (Table 3). The optimal resolution appeared to be at the lowest particle count (5 particles) for this specific prediction method, with the lowest mean error of 1.17 (see Table 3).

Particles	Prediction Error	Sample Size (N)
5	1.17	18
10	3.00	33
15	3.27	33
20	4.09	33
30	4.54	33
40	4.72	33
50	4.27	33

Discussion

Recapitulation of Study

This study demonstrates that dark matter subhaloes produce detectable gravitational perturbations in stellar streams, where gap counts show a strong positive correlation () with actual subhalo abundance across all tested stream resolutions. Several patterns can be spotted: a consistent baseline of false positive gaps in unperturbed streams with the specific false positive rate dependent on stream resolution, overprediction of subhalo density at low abundances (ratio 5.6 for single subhaloes) that approaches the accurate prediction (ratio 1) at higher subhalo counts, and the counterintuitive finding that prediction error increases with stream resolution.

Further, this study addresses all three primary research objectives: forward-modeling gravitational interactions across a defined parameter space, demonstrating that gap statistics can infer subhalo abundance (the inverse problem), and validating a model that produces physically consistent perturbations²⁸. The validation, however, also reveals crucial complexities, including a systematic prediction bias and a “Resolution Paradox,” which we now contextualize within the broader field.

Contextualization of Results

The strong positive correlation (r 0.811) between subhalo number and gap count confirms a fundamental relationship used to interpret perturbations in real streams (Figure 6). For example, the prominent “spur” and gap in the GD-1 stream were initially interpreted as potential evidence of a dark matter subhalo impact⁸. This study’s results provide further basis for this interpretation and demonstrates that such a signal can indeed arise from gravitational interactions alone. However, our observation of a consistent false-positive gap rate in unperturbed streams is a caveat to this correlation (Figure 4). This concern is substantiated by observational evidence: the complex S-shaped morphology of the Hyades tidal tails demonstrates that intricate structures can form through pure stellar dynamics in the Galactic tidal field, without any dark matter subhalo interactions²⁹. Consequently, to mitigate this false positive rate, any gap-counting method, including our heuristic, must contend with this inherent ambiguity when claiming a dark matter detection. Furthermore, our model’s systematic overprediction of subhalo density (especially at low subhalo abundances) highlights the major challenge of distinguishing dark matter-induced gaps from those caused by other effects (Table 2). Hence, there is an increasing importance of probabilistic approaches that assess the likelihood of gaps being caused by known baryonic perturbations (e.g., globular clusters) versus dark matter subhaloes³⁰.

The apparent contradiction between optimal correlation strength at 40 particles (Table 1) and maximal prediction error at the same resolution (Table 3) reflects the ‘Resolution Paradox.’ The strong correlation implies a high sensitivity for determining whether subhaloes are present, while the large prediction error reflects poor precision in quantifying exactly how many subhaloes are present; this is not a logical issue because the two metrics measure different tasks, but it is unexpected because one intuitively assumes that better data (in this case, higher resolution) should improve performance on both fronts simultaneously. This seeming “paradox” illuminates why the field is evolving beyond simple heuristics. In fact, modern methods are moving beyond simple gap-counting to account for the full complexity of stream morphology. For instance, a recent flexible framework that explicitly models off-track features like the GD-1 spur has been developed, fully acknowledging that a complete picture of the mechanical interactions requires more than just a density track³¹. Our “Resolution Paradox” mathematically demonstrates why this method advancement is necessary: as data quality improves (akin to our higher particle counts), a simple heuristic becomes inadequate, but the richer data simultaneously contains more complex signals that advanced methods can then exploit.

Finally, the subhalo-to-host mass ratio of 6% used in this simulation is a deliberate simplification to produce stronger and more interpretable signals in the stellar streams. This contrasts with the state-of-the-art observational efforts that aim to measure the Milky Way’s true mass distribution. The global model by Ibata et al. constraints the galaxy’s virial mass to —our work does not seek to replicate such precise values, but rather to validate the foundational physics principle of subhalo-gap correlation in a controlled environment³². The logical progression, as we suggest, is to build upon this foundation by incorporating a more cosmologically realistic subhalo mass function, thus enabling future versions of our model to contribute to quantitative estimates like those from the STREAMFINDER atlas.

Limitations

Taken altogether, these findings bring to light the potential promises and challenges of gap-based detection; however, they must also be interpreted with awareness of the numerous simplifying assumptions of this study. The uses of equal-mass subhaloes, a static galactic potential, and symmetric initial conditions do not fully depict the complexity of real Milky Way substructure. The omission of baryonic effects and relativistic dynamics, though reasonable for this proof-of-concept study, also limits direct comparison with observational data. Additionally, the relatively small parameter space (50 particles, 10 subhaloes) represents only a fraction of the expected substructure in galactic halos. Finally, as shown through the results, the gap-counting algorithm is useful for establishing baseline relationships but lacks the sophistication needed for more precise subhalo quantification—particularly for more complex environments.

More specifically, the gap-counting method is less sophisticated than two other emerging classes of methods. First, probabilistic gap-statistics models, such as those used to analyze the GD-1 stream³³, focus on calculating the likelihood that all observed gaps can be explained by known baryonic perturbers. Though our method infers subhalo density from gap counts, the probabilistic approach disfavors baryonic origins, and hence strengthens the case for dark matter interpretation of specific gaps. Second, a more recent work has used simulation-based inference (SBI) with machine learning, such as graph neural networks (GNNs), to perform field-level inference directly from a stream’s full 6-D phase-space information³⁴. These SBI methods do not rely on identifying gaps as an intermediate summary statistic; instead, they compress the entire complex morphology of the stream—including gaps, spurs, kinematic disturbances—to infer subhalo properties like mass and velocity. The resolution paradox we identified, where higher resolution improves correlation strength but also increases prediction error, is a direct consequence of our simple heuristic’s inability to decide complex interference patterns that these more holistic methods are designed to unveil.

Granted, our findings highlight the inherent limitations of gap-counting heuristics that point to a direction of future methodological developments. Rather than strictly refining simple gap-counting further, a more promising path forward involves adopting more sophisticated statistical frameworks that are emerging. Specifically:

1. Embrace Simulation-Based Inference (SBI): To mitigate the biases and errors of our heuristic method, future work should transition to SBI techniques which takes gap count into consideration, along with many other factors relating to a change in stream morphology.
2. Incorporate Probabilistic Modeling: Our model’s consistent overestimation of subhalo density at the low-density end reveals the limits of a deterministic gap-counting approach. A probabilistic framework could overcome this by combining prior knowledge from simulations with observed gap statistics, while also modeling the probability of false positives.
3. Expand Physical Realm: As a precursor or in conjunction with SBI, our model should be extended to include a realistic subhalo mass function and a time-dependent galactic potential. This is essential for producing simulations that are accurate enough for meaningful comparison to observational data from streams like GD-1 or Pal 5.

Conclusion

Our findings collectively reveal a fundamental tension in using stellar streams as dark matter detectors. Their long, extended structures and dynamical timescales make them highly sensitive to dark matter, but simultaneously, it is these same characteristics that create interference patterns that complicate seemingly straightforward gap-counting methods. Hence, future efforts may need to embrace this complexity and develop tools capable of interpreting the gravitational signals encoded in gaps and spirals to the fullest. Such approaches are essential to ultimately expose the long-hidden structure of the dark universe.

References

1. N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J.P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J.F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J.M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, F. K. Hansen, D. Herranz, S. R. Hildebrandt, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J.M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis , M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y.Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M.A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J.L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A.S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, A. Zonca. Planck 2018 results. VI. Cosmological parameters. A&A. 641, A6 (2020). [↩]
2. W. Jiang, J. Han, F. Dong, F. He. Self-similar decomposition of the hierarchical merger tree of dark matter halos. The Astrophysical Journal. 988, 160. (2025). [↩]
3. J. J. Webb, J. Bovy. And In The Darkness Unbind Them: High-Resolution Simulations of Dark Matter Subhalo Disruption in a Milky Way-like Tidal Field. Monthly Notices of the Royal Astronomical Society. 499, 116–128 (2020). [↩]
4. R. Massey, T. Kitching, J. Richard. The dark matter of gravitational lensing. Reports on Progress in Physics. 73, 086901 (2010). [↩]
5. N. Dalal, Y.D. Hezaveh, D.P. Marrone, Y.Y. Mao, W. Morningstar, D. Wen, R.D. Blandford, J.E. Carlstrom, C.D. Fassnacht, G.P. Holder, A. Kemball, P.J. Marshall, N. Murray, L. Perreault Levasseur, J.D. Vieira, R.H. Wechsler. DETECTION OF LENSING SUBSTRUCTURE USING ALMA OBSERVATIONS OF THE DUSTY GALAXY SDP.81. ApJ. 823, 37 (2016). [↩]
6. P. Menker, A. Benson. Advancing Stellar Streams as a Dark Matter Probe — I: Evolution of the CDM subhalo population. arXiv:2406.11989 (2024). [↩]
7. A. Bonaca, A. M. Price-Whelan. Stellar Streams in the Gaia Era. New Astronomy Reviews. 100, 101713 (2024). [↩] [↩]
8. A. Bonaca, D.W. Hogg, A.M. Price-Whelan, C. Conroy. The Spur and the Gap in GD-1: Dynamical evidence for a dark substructure in the Milky Way halo. ApJ. 880, 38 (2019). [↩] [↩] [↩]
9. D. Erkal, V. Belokurov. Forensics of subhalo–stream encounters: the three phases of gap growth. Monthly Notices of the Royal Astronomical Society. 450, 1136–1149 (2015). [↩]
10. J. Sanders, J. Bovy, D. Erkal. Dynamics of stream-subhalo interactions. Monthly Notices of the Royal Astronomical Society. 457, 3817–3835 (2015). [↩]
11. V. Springel, S. D. M. White, A. Jenkins, C.S. Frenk, N. Yoshida, L. Gao, J. Navarro, R. Thacker, D. Croton, J. Helly, J.A. Peacock, S. Cole, P. Thomas, H. Couchman, A. Evrard, J. Colberg, F. Pearce. Simulations of the formation, evolution and clustering of galaxies and quasars. Nature. 435, 629–636 (2005). [↩]
12. J. Hermans, N. Banik, C. Weniger, G. Bertone, G. Louppe. Towards constraining warm dark matter with stellar streams through neural simulation-based inference. Monthly Notices of the Royal Astronomical Society. 507, 1999–2011 (2021). [↩]
13. S.R. Green, R.M. Wald. Newtonian and Relativistic Cosmologies. Phys. Rev. D. 85, 063512 (2012). [↩]
14. V. Springel, S. D. M. White, A. Jenkins, C.S. Frenk, N. Yoshida, L. Gao, J. Navarro, R. Thacker, D. Croton, J. Helly, J.A. Peacock, S. Cole, P. Thomas, H. Couchman, A. Evrard, J. Colberg, F. Pearce. Simulations of the formation, evolution and clustering of galaxies and quasars. Nature. 435, 629–636 (2005). [↩]
15. A. Eyre, J. Binney. The mechanics of tidal streams. Monthly Notices of the Royal Astronomical Society. 413, 1852–1874 (2011). [↩]
16. E.O. Nadler, S. Birrer, D. Gilman, R.H. Wechsler, X. Du, A. Benson, A.M. Nierenberg, T. Treu. Dark Matter Constraints from a Unified Analysis of Strong Gravitational Lenses and Milky Way Satellite Galaxies. ApJ. 917, 7 (2021). [↩]
17. R.W Robinett. Using Physics to Learn Mathematica to Do Physics: From Homework Problems to Research Examples. arXiv:0712.2358 (2007). [↩]
18. T. Quinn, N. Katz, J. Stadel, G. Lake. Time stepping N-body simulations. arXiv:astro-ph/9710043 (1997). [↩]
19. W. Dehnen. A Hierarchical O(N) Force Calculation Algorithm. Journal of Computational Physics. 179, 27–42 (2002). [↩]
20. J. E. Barnes. Gravitational softening as a smoothing operation. Monthly Notices of the Royal Astronomical Society. 425, 1104–1120 (2012). [↩]
21. A. Shirokov. Gravitational Softening and Adaptive Mass Resolution. arXiv:0711.2989 (2008). [↩]
22. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press (2007). [↩]
23. J. Eckmann, F. Hassani. The detection of relativistic corrections in cosmological N-body simulations. Celest. Mech. Dyn. Astron. 132, 1-2 (2020). [↩]
24. J. Bovy, D. Erkal, J.L. Sanders. Linear perturbation theory for tidal streams and the small-scale CDM power spectrum. Monthly Notices of the Royal Astronomical Society. 466, 628–668 (2017). [↩] [↩]
25. E. Contini, G.D. Lucia, S. Borgani. Statistics of substructures in dark matter haloes. Monthly Notices of the Royal Astronomical Society. 420, 2978-2989 (2012). [↩]
26. A. Bonaca, A. M. Price-Whelan. Stellar Streams in the Gaia Era. New Astronomy Reviews. 100, 101713 (2024). [↩]
27. D. Erkal, V. Belokurov, J. Bovy, J.L. Sanders. The number and size of subhalo-induced gaps in stellar streams. Monthly Notices of the Royal Astronomical Society. 463, 102–119 (2016). [↩]
28. J. Lu, T. Lin, M. Sholapurkar, A. Bonaca. Detectability of dark matter subhalo impacts in Milky Way stellar streams. arXiv:2502.07781 (2025). [↩]
29. S. Meingast, J. Alves. Extended stellar streams in the solar neighborhood. A&A. 621, 6 (2019). [↩]
30. S. Ferrone, M. Montuori, P.D. Matteo, A. Mastrobuono-Battisti, R. Ibata, P. Bianchini, S. Khoperskov, N. Leclerc, C. Hottier, E. Stein, D. Valls-Gabaud, O.N. Snaith, M. Haywood. Gaps in stellar streams as a result of globular cluster fly-bys. A&A. 699, A289 (2025). [↩]
31. K. Tavangar, A.M. Price-Whelan. Inferring the density and membership of stellar streams with flexible models: The GD-1 stream in Gaia Data Release 3. The Astrophysics Journal. 988, 45 (2025). [↩]
32. R. Ibata, K. Malhan, W. Tenachi, A. Ardern-Arentsen, M. Bellazzini, P. Bianchini, P. Bonifacio, E. Caffau, F. Diakogiannis, R. Errani, B. Famaey, S. Ferrone, N. Martin, P. di Matteo, G. Monari, F. Renaud, E. Starkenburg, G. Thomas, A. Viswanathan, Z. Yuan. Charting the Galactic acceleration field II. A global mass model of the Milky Way from the STREAMFINDER Atlas of Stellar Streams detected in Gaia DR3. The Astrophysics Journal. 967, 89 (2024). [↩]
33. Y. Doke, K. Hattoru. Probability of forming gaps in the GD-1 stream by close encounters of globular clusters. ApJ. 941, 129 (2022). [↩]
34. P.X. Ma, K.K. Rogers, T.S. Li, R. Hložek, J.J. Webb, R. Huang, J.Meunier. Towards characterizing dark matter subhalo perturbations in stellar streams with graph neural networks. ApJ. 987, 96 (2025). [↩]