Log In Sign Up

Sensitivity Estimation for Dark Matter Subhalos in Synthetic Gaia DR2 using Deep Learning

by   Abdullah Bazarov, et al.

The abundance of dark matter subhalos orbiting a host galaxy is a generic prediction of the cosmological framework. It is a promising way to constrain the nature of dark matter. Here we describe the challenges of detecting stars whose phase-space distribution may be perturbed by the passage of dark matter subhalos using a machine learning approach. The training data are three Milky Way-like galaxies and nine synthetic Gaia DR2 surveys derived from these. We first quantify the magnitude of the perturbations in the simulated galaxies using an anomaly detection algorithm. We also estimate the feasibility of this approach in the Gaia DR2-like catalogues by comparing the anomaly detection based approach with a supervised classification. We find that a classification algorithm optimized on about half a billion synthetic star observables exhibits mild but nonzero sensitivity. This classification-based approach is not sufficiently sensitive to pinpoint the exact locations of subhalos in the simulation, as would be expected from the very limited number of subhalos in the detectable region. The enormous size of the Gaia dataset motivates the further development of scalable and accurate computational methods that could be used to select potential regions of interest for dark matter searches to ultimately constrain the Milky Way's subhalo mass function.


Cataloging Accreted Stars within Gaia DR2 using Deep Learning

The goal of this paper is to develop a machine learning based approach t...

Simulation-Assisted Decorrelation for Resonant Anomaly Detection

A growing number of weak- and unsupervised machine learning approaches t...

Modeling Representation of Videos for Anomaly Detection using Deep Learning: A Review

This review article surveys the current progresses made toward video-bas...

Detecting Attacks on IoT Devices using Featureless 1D-CNN

The generalization of deep learning has helped us, in the past, address ...

Simulation Assisted Likelihood-free Anomaly Detection

Given the lack of evidence for new particle discoveries at the Large Had...

Inferring dark matter substructure with astrometric lensing beyond the power spectrum

Astrometry – the precise measurement of positions and motions of celesti...

From Dark Matter to Galaxies with Convolutional Networks

Cosmological surveys aim at answering fundamental questions about our Un...

1 Introduction

Dark matter (DM) represents roughly 84% of the matter content in the Universe (Planck Collaboration et al., 2020). However, unveiling its nature has proven a difficult endeavour, and none of the proposed candidates (from several extensions of the Standard Model to primordial black holes) have yet been detected. Cold DM is expected to form subhalos with masses many orders of magnitude below (Blumenthal et al., 1984; Wang et al., 2020a), which is roughly the mass above which galaxies can form (Zavala and Frenk, 2019). The abundance of subhalos is dependent on the nature of DM. This dependency can be explained, on the one hand, by the effect of the properties of the DM on the linear matter power spectrum. If for cold DM the minimum halo mass might be as small as (Zybin et al., 1999; Bringmann, 2009), microscopic properties of the DM particle, e.g. non-negligible thermal velocities or quantum pressure, introduce a cut-off at the small scales in alternative DM scenarios. On the other hand, the nature of DM, e.g. self-interactions, further impacts the non-linear growth of structures (Schneider et al., 2012; Vogelsberger et al., 2016). Detecting a dark subhalo would be the first direct evidence of DM clustering at small scales. Furthermore, constraints on the subhalo abundance would provide valuable information about the particle nature of DM.

Subhalos with masses lower than are unable to form stars and remain dark, thus hindering their detection. Strategies that aim to detect dark subhalos rely on measuring their gravitational signatures via stellar dynamics (Ibata et al., 2002; Yoon et al., 2011; Carlberg, 2012; Bovy et al., 2017; Banik et al., 2018; Bonaca et al., 2019; Benito et al., 2020; Feldmann and Spolyar, 2015; Buschmann et al., 2018) or gravitational lensing (Hezaveh et al., 2016; Van Tilburg et al., 2018; Díaz Rivero et al., 2018; Gilman et al., 2019; Brehmer et al., 2019; Vattis et al., 2020) and, in the case of several DM candidates, e.g. Weakly Interacting Massive Particles (WIMPs), on detecting the flux of final stable particles produced by DM annihilation or decay (Moliné et al., 2017; Calore et al., 2019; Coronado-Blázquez et al., 2019, 2021; Mirabal and Bonaca, 2021). The goal of searches based on stellar dynamics is to detect perturbations in the phase-space distribution of Milky Way (MW) stars induced by gravitational effects of passing subhalos. We can look for these perturbations in stellar streams (Ibata et al., 2002; Yoon et al., 2011; Carlberg, 2012; Bovy et al., 2017; Banik et al., 2018; Bonaca et al., 2019) and in the disk or the halo stars (Feldmann and Spolyar, 2015; Buschmann et al., 2018). In the present work we investigate the usage of an anomaly detection and classification algorithms in the search for the imprint caused in halo stars by passing substructures. In this way, we exploit the increasing size of observational datasets and state-of-the-art techniques in deep learning.

In recent years, deep learning techniques have been applied in the search for substructures in our Galaxy (Ostdiek et al., 2020; Necib et al., 2020; Shih et al., 2021). These detection methods assume that stars in the MW sharing a common origin should cluster in orbital properties and/or composition. Our search differs in that we aim to identify stars that, regardless of their origin, have their distribution in phase-space perturbed by the passage of a dark matter subhalo. For any identified star, it must be possible to test the halo hypothesis independently of the methodology used to select the candidates. One possibility could be to preselect the stars using a ML

-based classifier, followed by detailed hypothesis tests using e.g. the orbital arc method 

(Kipper et al., 2020, 2021) or the stellar wakes technique (Buschmann et al., 2018).

The raw data that we used, which are described in section 2, are three MW-like galaxies from the Latte suite of FIRE-2 simulations (Wetzel et al., 2016) and nine synthetic Gaia DR2 surveys generated from the simulated galaxies by means of the Ananke framework (Sanderson et al., 2020). First, we processed the synthetic Gaia datasets to correlate the position of stars and the dark subhalos, which were previously identified in the simulated galaxies. In section 3, we estimate the detectability of the subhalo-associated stars using deep learning techniques. We conclude in section 4.

2 Datasets

As our raw data, we have used three MW-like galaxies from the Latte suite of FIRE-2 simulations (Wetzel et al., 2016; Garrison-Kimmel et al., 2017; Hopkins et al., 2018) (dubbed m12f, m12i and m12m) and nine synthetic Gaia DR2 surveys (Sanderson et al., 2020). This section describes these datasets and the processing we have performed on them.

2.1 Milky Way-like Galaxies

We have used the simulation snapshots at of three MW-like galaxies, namely m12f, m12i and m12m (Wetzel et al., 2016; Garrison-Kimmel et al., 2017; Hopkins et al., 2018). In the following we briefly describe how these MW analogues were obtained. For a complete description of this and the details of the N-body simulations we refer the interested reader to (Wetzel et al., 2016) and references therein. The MW analogues were first identified in a DM-only cosmological simulation requiring that at : (i) their virial mass is in the range of 111Virial mass and virial radius follow the relation , with the average matter density of the Universe. (which agrees with recent measurements Wang et al. 2020b; Karukes et al. 2020; Shen et al. 2021) and (ii) there is no neighboring halo of similar mass within . Three halos selected in this manner were then simulated using the zoom-in technique (Oñorbe et al., 2014). Simulations were run using the Gizmo gravity plus hydrodynamics code in meshless finite-mass (MFM) mode (Hopkins, 2015) and the FIRE-2 baryonic physics model (Hopkins et al., 2018).

We have identified DM subhalos in snapshots at of the MW-like galaxies using the Amiga Halo Finder (AHF) code (Knollmann and Knebe, 2009). The AHF algorithm identifies bound DM

structures by hierarchically clustering 3D positions of

DM particles in the simulation. Following (Garrison-Kimmel et al., 2017), AHF was run only on DM particles. We have selected subhalos with more than 85 DM particles (corresponding to subhalos with masses ) since those substructures are reliably resolved in the simulation (Garrison-Kimmel et al., 2017). Approximately subhalos222AHF identifies 1298, 1001 and 1281 subhalos with for m12f, m12i and m12m, respectively. for each MW-like galaxy remain as potentially observable. Figure 1 shows the cumulative subhalo mass function normalized by the virial mass of the host halo (left panel) and the radial distribution of the subhalo population normalized by the virial radius (right panel). The virial masses are , and for m12f, m12i and m12m, respectively.333In our work we have used the virial mass definition given by , with the critical density of the Universe at . In figure 2 we show the mass of the subhalos as a function of their galactocentric distance for m12f, m12i and m12m. It should be noted that no subhalos are identified below from the center of the galaxies, as previously noted in (Garrison-Kimmel et al., 2017). Furthermore, the 97%, 91% and 94% of the subhalos below 50 kpc for m12f, m12i and m12m, respectively, have masses lower than . The most massive subhalo below 50 kpc is identified at 20 kpc with for m12f, at 43 kpc with for m12i and at 43 kpc with for m12m. The depletion of the most massive dark subhalos in the inner 50 kpc of the MW-like galaxies conditions the ability to identify stars in the stellar halo, which have been perturbed by the passage of a dark matter subhalo, using the deep-learning techniques explored in this study.

Figure 1: The subhalo mass function scaled to the host halo virial mass (left) and the radial distribution of the subhalo population (right) for the three MW-like galaxies used in this work.
Figure 2: Subhalo mass as a function of galactocentric distance. Each dot corresponds to a subhalo as identified by AHF for m12f (left), m12i (middle) and m12m (right) galaxies.

2.2 Synthetic Gaia Surveys

The nine synthetic Gaia DR2 surveys were generated by applying the Ananke framework (Sanderson et al., 2020) to the three MW-like galaxies. Per simulated galaxy, three synthetic surveys were generated by adopting three local standards of rest. Each synthetic survey contains approximately a billion mock stellar observations resembling Gaia DR2. We restrict our attention to stellar halo stars, applying a selection in true vertical distances kpc. In this way, we remove disk stars that could suffer from disturbances induced, for example, by spiral arms, the Galactic bar or giant molecular clouds. We are left with mock stars for each LSR for the subsequent analysis. This reduced dataset consists of nearly 2 billion observed stars for the three different MW-like galaxies, three LSRs for each, correlated with potentially observable DM subhalo locations. Table 1 summarizes some statistics of this dataset.

Figure 3: Number of stars associated to a particular subhalo as a function of its mass for m12f (left), m12i (middle) and m12m (right) galaxies.

Stars are tagged as halo-associated if their true distance to the central position of a subhalo is lower than 1 kpc. It is to be noticed that these halo-associated stars might not be bound to the subhalos. Figure 3 shows the total number of stars associated to a subhalo as a function of the subhalo’s mass for each LSR and each simulated galaxy. Within each galaxy, less than of the subhalos contain associated stars, and 66%, 84% and 75% of this fraction contain less than 10 associated stars for m12f, m12i and m12m galaxies, respectively. Furthermore, approximately 40% of the halos that have associated stars contain only one star. The m12f galaxy has a larger percentage of subhalos which are associated to more than 100 stars compared to that of m12i or m12f galaxies. This is because the former galaxy has a larger fraction of subhalos below 30 kpc. We plot the projected stellar number densities for all LSRs on figure 4, along with the halo locations and halo-associated observed stars.

stars with kpc halo-associated stars [] with [] subhalos w/ halo-associated stars
m12f LSR0 216,446,024 0.0291% 0.35% 73
LSR1 182,538,592 0.0291% 0.32% 76
LSR2 204,017,261 0.0306% 0.35% 71
m12i LSR0 139,167,343 0.0019% 0.41% 63
LSR1 132,655,442 0.0017% 0.41% 61
LSR2 131,474,668 0.0010% 0.23% 67
m12m LSR0 170,255,144 0.0013% 0.09% 67
LSR1 156,093,757 0.0016% 0.12% 71
LSR2 161,369,511 0.0013% 0.19% 68
Table 1: Summary statistics of synthetic Gaia DR2 reduced catalogs used in this work (see text for more details).
Figure 4: Projected stellar number densities in the synthetic Gaia datasets in true galactocentric coordinates for the three MW-like galaxies, namely m12f (top row), m12i (middle row) and m12m (bottom row) for LSR0. The halo locations are shown in black, with the size of the markers being proportional to the halo’s virial radius, while the regions with halo-associated observed stars are shown in red. The Sun’s position for each case is .

3 Deep Learning Search of Subhalo-associated Stars

Dark subhalos perturb the positions and velocities of nearby stars. We wish to estimate if these imprints are detectable in MW-like galaxies and in synthetic data that accounts for observational uncertainties. Let us assume, without loss of generality, that the properties

of each star particle (or observed star) are drawn from the probability distribution


if the star particle (observed star) has or has not been affected, respectively, by a dark subhalo at a given time in the Latte (Ananke) simulation. Then, if the probabilities are known, the likelihood ratio

is the optimal discriminator between the two hypotheses for a given observation according to the Neyman-Pearson lemma (Neyman and Pearson, 1992). These probabilities are not known, however, and we only have simulated examples of either halo-associated or background star particles (observed stars). In the following, we investigate the possibility of using machine learning to define an approximate discriminator between the two hypotheses, and thus quantify the difference between the halo-associated stars and the background.

3.1 Detectability

As a starting point, we first focus on the Latte simulations, where for each star particle, the full six-dimensional phase-space coordinates, namely the three-dimensional Galactocentric Cartesian positions and velocities, are known. Unlike for the synthetic Gaia dataset, the disk is not excluded at this stage. For each star particle, we compute the Euclidean distance to the nearest dark subhalo , and if it is below a threshold  kpc, we identify the star particle as a halo-associated or signal particle. We then use an anomaly detection approach to estimate the strength of the subhalo signal (Baldi and Hornik, 1989; Sakurada and Yairi, 2014). For this purpose, the background-only likelihood

is indirectly approximated using a so-called autoencoder neural network, and deviations from the background-only distribution are quantified.

Each star particle is characterized by the feature vector

containing its three-dimensional position and velocity, i.e. . Let us define an encoder and a decoder as


respectively, such that approximates for any given input via a lower-dimensional representation. Both the encoder and decoder are implemented as feedforward neural networks, optimized by tuning the weights using only the background examples as follows:


By construction, the encoder-decoder will tend to reconstruct well the background-like samples that it was optimized on. On the other hand, for any other that is not distributed as , we would expect on average higher values for the reconstruction loss . Therefore, we can use the distribution , optimized only on the background particles, as an empirical discriminator between the background and signal samples. We have checked this approach by defining a fake signal consisting of a random sub-population of stars irrespective of the dark subhalo locations. In this case, no detectable difference between the main sample and the random subpopulation of stars is expected with this method.

We optimize the model on the m12m galaxy, while cross-checking the performance on the m12i

galaxy. This ensures that the model is not simply memorizing the locations of the halos, as the result in this case would not be generalizable to other galaxy simulations. Training is carried out for 100 iterations (epochs) over the full dataset using the Adam optimizer 

(Kingma and Ba, 2014) with a learning rate of and a minibatch size of star particles. We show the evolution of the total reconstruction loss over training epochs on figure 5 for both the real subhalo signal (left panel) and the fake signal cases (right panel).

Figure 5: Relative reconstruction loss with respect to the first epoch for the training (m12m) and testing (m12i) datasets based on the anomaly detection model over the training epochs. By construction, we expect the loss in m12m to decrease as the model fits the samples. The training is performed on star particles excluding the subhalo-associated particles (left) and by excluding the same number of random particles (right) as a check.

We observe that the model converges for m12m and exhibits in general stable behaviour for m12i.

Figure 6 shows the distributions for the signal and background stars for the m12f dataset never used for training. It is observed that for the real subhalo signal case (left panel), there is a distinction between the distribution of halo-associated (signal) and non-halo associated (background) stars, with the signal stars having on average higher values of the reconstruced distribution. No such distinction is observed for the model trained and tested on the random subset (right panel), as would be expected.

Figure 6: The distribution of the reconstruction loss for the m12f galaxy that was never used in the training procedure. We show the distributions of the halo-associated (signal) and the rest (background) of the star particles separately. On the left, the model was trained and evaluated with the subhalo-associated as the signal, while on the right, a random set of stars was denoted as the signal as a check.

We quantify the performance of the anomaly detection in terms of the true positive and false negative rates. The true positive rate (TPR) gives the fraction of signal stars that are correctly identified as signal particles at a particular threshold (i.e. a given value of ),


Contrary, the false positive rate (FPR) is the fraction of background stars that are incorrectly identified as signal, namely


Figure 7 shows the FPR versus TPR while scanning over for the real (solid blue) and fake (dashed black) signal cases.

Figure 7: The true positive rate (horizontal axis) vs. the false positive rate (vertical axis) when the reconstruction loss is used as a discriminator between the signal and background labels on the star particles from the m12f galaxy. By construction, we observe no distinction for random stars, while the model distinguishes halo-associated stars from the background using a combination of the position and velocity information.

For the real or subhalo case, we see that at a , the FPR is (i.e. 80% of the signal stars are correctly identified while we misclassify 15% of the background stars as signal), presenting a significant improvement over a random selection.

Based on this statistical anomaly detection model optimized on the m12m simulated galaxy, we can conclude that halo-associated star particles in the m12f dataset have a distinguishable distribution in 6-dimensional phase-space consisting of positions and momenta and that the anomaly detection method is able to correctly identify halo-associated stars.

3.2 Feasibility in Synthetic Gaia Survey

In this section, we investigate if the subhalo-associated stars are detectable in the synthetic Gaia survey derived from this simulation, that is, under the effects of extinction, partial measurement of the radial velocity and measurement errors. This is done by searching for dark subhalos on the reduced synthetic surveys described in section 2.2. The goal is again to select candidate stars which are likely to be perturbed by a nearby dark subhalo, such that they could potentially be further analyzed with more detailed approaches. The search for dark subhalos in the reduced Gaia-like catalogs differs from the previous section on two fronts. On the one hand, the mock observed stars in the synthetic Gaia datasets are divided into patches using the hierarchical pixelization algorithm HEALPY (Zonca et al., 2019; Górski et al., 2005) with a pixel level 6. This allows to process the data in manageable subsets in a physically meaningful fashion. In addition, as the subhalo-associated stars are located in well-defined, localized regions in the sky, we avoid using the absolute right ascension and declination coordinates to unfairly bias the model. Instead, we compute the positional information with respect to the pixel center. For a Gaia DR2-like dataset, a pixel can contain up to stars.

On the other hand, the input feature of each observed star is different. For each synthetic dataset realization , , we then have a list of (star observation, label) pairs

Each stellar feature vector consists of the following astrometric observables:

  • parallax [mas],

  • the right ascension with respect to the pixel center [deg],

  • the declination with respect to the pixel center [deg],

  • the proper motion in the right ascension direction (multiplied by ) [mas/year],

  • the proper motion in the declination [mas/year] and

  • the radial velocity [km/s].

These observables are available with estimated uncertainties, resulting in 12 input features. Features which are not always measured, such as the radial velocities, are filled with placeholder values for a consistent numerical treatment.

For the synthetic Gaia dataset, we cross-check the anomaly detection approach described in the previous section against a simple binary classifier, where the signal model is used explicitly, but which is thus limited by the available statistics for the halo-associated stars. Contrary, in the anomaly-detection based approach, the star labels based on the proximity to a dark subhalo were only used to exclude the signal samples from the optimization. The supervised classification model uses the signal sample labels directly, i.e. the optimization target is the star label for background and signal stars, respectively. Thus, it can be used to determine the upper limit detectability for this particular signal model, assuming training statistics are not a limiting factor.

The binary star classification model is defined as a parametric function using a deep neural network


which can be optimized by tuning the weights

to minimize a classification loss function. We use the the focal loss, which is a modification of the binary cross entropy loss, originally proposed for rare object detection in 

(Lin et al., 2017), and is defined as the following sum over the total number of pixels in the dataset:


where and are empirical factors that adjust the weight of easy-to-classify background-like examples in the loss. We choose based on the defaults introduced in (Lin et al., 2017). By this construction, the model output for star

is a continuous value between 0 and 1 that can be interpreted as a test statistic for the star being labeled as signal. We use dropout regularization 

(Srivastava et al., 2014) with a probability of in the training phase to limit the amount of overtraining. The anomaly detection model was trained for 200 epochs, while the classification model for 50 epochs. As before, we use m12m for the optimization, m12i for testing, while m12f is used for the final physics validation. The training and testing is done on stars from all three LSRs simultaneously. Overall, as summarized in Table 1, the optimization, testing and final validation is carried out on nearly 1.5 billion mock stars, of which less than 0.01% are identified as signal, resulting in extreme class imbalance as well as an overall low number of independent signal samples.

The sensitivity of the anomaly detection and classifier methods for identifying halo-associated stars in the synthetic Gaia dataset can be seen on figure 8. As before, the m12f dataset was never used in the optimization. We observe that the binary classification distinguishes between the halo-associated and background stars at a non-negligible level, with a FPR of at a TPR of . On the other hand, the anomaly detection approach, where we only attempt to learn the background distribution, does not differ significantly from a purely random selection in the synthetic survey.

Figure 8: The true positive rate versus the false positive rate after evaluating the anomaly-detection (blue) and binary classification (orange) models on the synthetic Gaia dataset m12f.

Based on this, we conclude that the detectability of the halo-associated stars from the astrometric properties in the synthetic Gaia-DR2 surveys is not straightforward, even if the dark matter subhalos have a detectable kinematic effect on the underlying simulated star particles, as demonstrated in the previous subsection.

4 Conclusions

ML techniques, either alone or combined with classical methods, have been demonstrated to be helpful in uncovering new structures in Gaia-scale datasets (e.g. Necib et al. 2020). Dark subhalos are among the most challenging substructures to search for. In this paper, we study the detectability of dark subhalos by means of ML in three MW-like galaxies and in nine synthetic Gaia DR2 surveys. Rather than attempting to pinpoint the exact subhalo locations and determine their properties, we attempt to identify candidate stars that are likely to be close to a subhalo on a statistical basis.

We have first correlated star particles in the simulated galaxies and mock stars in the synthetic catalogs with the position of dark subhalos found by the AHF. In section 3.2 we then tested against simulated galaxies the feasibility of an anomaly detection algorithm to detect the phase-space imprint in stellar halo stars of nearby subhalos. This algorithm builds a likelihood function of the background star particles and is able to correctly identify 80% of signal stars while misclassifying as signal 15% of background particles. We concluded that the distribution function of the 6-dimensional phase-space coordinates of signal and background star particles are distinguishable in the MW-like galaxies used in this work. Therefore, on a statistical basis, the position and velocity information can be combined into a statistical discriminator for the halo-associated signal.

Finally, we have tested the feasibility of the anomaly detection algorithm in Gaia DR2-like surveys in section 3.2. The unsupervised anomaly detection algorithm was cross-checked against a binary supervised classification algorithm that uses the signal model explicitly. The anomaly detection approach has no sensitivity to distinguish between signal and background stars, while the binary classification algorithm is able to select 50% of signal stars while wrongly identifying 15% of background stars as signal. Although the binary classification shows a mild sensitivity, overall both approaches are of limited effectiveness in the synthetic Gaia survey. The two main factors limiting the sensitivity of deep learning searches for stellar perturbations in phase-space induced by the passage of a DM subhalo are

  1. the low signal-to-noise ratio. Although the synthetic catalogs contain more stars than the real Gaia-DR2 

    (Sanderson et al., 2020), a very limited number of subhalos are contained in the region observed by the survey. Thus resulting in only a handful of halos being usable for the ML optimization. Furthermore, for practicality, we so far focused on the halo stars, neglecting the galactic disk, which is the bulk of the survey.

  2. There is not enough precision in mock observations. On top of the small number (according to the subhalo populations found for the MW-like galaxies in this analysis) of dark matter subhalos within 40 kpc from the Galactic center (GC), those subhalos have a mass below . The velocity change experienced by a star due to the encounter with a subhalo of is roughly , where v is the star’s initial velocity (Feldmann and Spolyar, 2015). Therefore, in order to measure kinematic perturbations the velocity precision must be . This level of precision might be achieved by the 4MOST Milky Way Halo Low-Resolution survey (Helmi et al., 2019).

Another possible limiting factor for deep searches of DM subhalos is the derivation of the synthetic surveys used for training. From each star particle in the Latte simulation, a set of observed stars are derived by sampling from its six-dimensional phase-space distribution (Sanderson et al., 2020). The sampling process introduces a smearing scale which distorts and possibly diminishes the imprint of the passage of a dark subhalo. We checked for this effect by applying the supervised classification approach to two small subsets: one only containing the observed stars whose true position and velocity equal that of the parent star, and another subset with the same fraction of stars that were randomly selected from those that have different true position and velocity with respect to the parent star particle from which they have been generated. In both cases, we obtain that the supervised classification algorithm has no sensitivity to distinguish between background and signal stars and the results are qualitatively the same for both subsets showing that we are limited by the available statistics. Any further study of the impact of the smearing process is left to future investigations.

A number of subsequent improvements to the methodology are possible. In the above analysis, all the observed stars were treated independently of each other. Local correlations, density or clustering were not taken into account, which could potentially limit the sensitivity of the method used so far. As an example, novel approaches based on density-based clustering have been employed for open clusters (Castro-Ginard et al., 2018) and may be interesting to study here for dark subhalos. Clustering can also be combined with unsupervised deep learning for anomaly detection (Mikuni and Canelli, 2021). Another possible approach is the direct search for overdensities by comparing signal and sideband regions, which has been so far demonstrated for stellar streams, but could potentially be studied also for dark subhalos (Shih et al., 2021). Furthermore, in order to understand the potential sensitivity of the method, simulated datasets with a known DM distribution were used. However, the halo distribution in these is fixed, and the number of actual simulated halos in the potentially visible region is limited. Additional simulated datasets with a varying halo distribution could be helpful to establish sensitivity dependence of a potential method on halo mass and distance from the galactic center.


This work was supported by the Estonian Research Council grants PSG700, PRG803, PRG1006, PUTJD907 and MOBTT5, MOBTP187, and by the European Regional Development Fund through the CoE program grant TK133.


  • Baldi and Hornik (1989) Baldi, P., Hornik, K., 1989.

    Neural networks and principal component analysis: Learning from examples without local minima.

    Neural Netw. 2, 53–58. URL:, doi:10.1016/0893-6080(89)90014-2.
  • Banik et al. (2018) Banik, N., Bertone, G., Bovy, J., Bozorgnia, N., 2018. Probing the nature of dark matter particles with stellar streams. JCAP 2018, 061. doi:10.1088/1475-7516/2018/07/061, arXiv:1804.04384.
  • Benito et al. (2020) Benito, M., Criado, J.C., Hütsi, G., Raidal, M., Veermäe, H., 2020. Implications of Milky Way substructures for the nature of dark matter. Phys. Rev. D 101, 103023. doi:10.1103/PhysRevD.101.103023, arXiv:2001.11013.
  • Blumenthal et al. (1984) Blumenthal, G.R., Faber, S.M., Primack, J.R., Rees, M.J., 1984. Formation of galaxies and large-scale structure with cold dark matter. Nature 311, 517--525. doi:10.1038/311517a0.
  • Bonaca et al. (2019) Bonaca, A., Hogg, D.W., Price-Whelan, A.M., Conroy, C., 2019. The Spur and the Gap in GD-1: Dynamical Evidence for a Dark Substructure in the Milky Way Halo. The Astrophysical Journal 880, 38. doi:10.3847/1538-4357/ab2873, arXiv:1811.03631.
  • Bovy et al. (2017) Bovy, J., Erkal, D., Sanders, J.L., 2017. Linear perturbation theory for tidal streams and the small-scale CDM power spectrum. MNRAS 466, 628--668. doi:10.1093/mnras/stw3067, arXiv:1606.03470.
  • Brehmer et al. (2019) Brehmer, J., Mishra-Sharma, S., Hermans, J., Louppe, G., Cranmer, K., 2019. Mining for Dark Matter Substructure: Inferring Subhalo Population Properties from Strong Lenses with Machine Learning. The Astrophysical Journal 886, 49. doi:10.3847/1538-4357/ab4c41, arXiv:1909.02005.
  • Bringmann (2009) Bringmann, T., 2009. Particle models and the small-scale structure of dark matter. New Journal of Physics 11, 105027. doi:10.1088/1367-2630/11/10/105027, arXiv:0903.0189.
  • Buschmann et al. (2018) Buschmann, M., Kopp, J., Safdi, B.R., Wu, C.L., 2018. Stellar Wakes from Dark Matter Subhalos. Phys. Rev. Letters 120, 211101. doi:10.1103/PhysRevLett.120.211101, arXiv:1711.03554.
  • Calore et al. (2019) Calore, F., Hütten, M., Stref, M., 2019. Gamma-Ray Sensitivity to Dark Matter Subhalo Modelling at High Latitudes. Galaxies 7, 90. doi:10.3390/galaxies7040090, arXiv:1910.13722.
  • Carlberg (2012) Carlberg, R.G., 2012. Dark Matter Sub-halo Counts via Star Stream Crossings. The Astrophysical Journal 748, 20. doi:10.1088/0004-637X/748/1/20, arXiv:1109.6022.
  • Castro-Ginard et al. (2018) Castro-Ginard, A., Jordi, C., Luri, X., Julbe, F., Morvan, M., Balaguer-Núñez, L., Cantat-Gaudin, T., 2018. A new method for unveiling open clusters in gaia-new nearby open clusters confirmed by dr2. Astronomy & Astrophysics 618, A59.
  • Coronado-Blázquez et al. (2021) Coronado-Blázquez, J., Doro, M., Sánchez-Conde, M.A., Aguirre-Santaella, A., 2021. Sensitivity of the Cherenkov Telescope Array to dark subhalos. Physics of the Dark Universe 32, 100845. doi:10.1016/j.dark.2021.100845, arXiv:2101.10003.
  • Coronado-Blázquez et al. (2019) Coronado-Blázquez, J., Sánchez-Conde, M.A., Di Mauro, M., Aguirre-Santaella, A., Ciucă, I., Domínguez, A., Kawata, D., Mirabal, N., 2019. Spectral and spatial analysis of the dark matter subhalo candidates among Fermi Large Area Telescope unidentified sources. JCAP 2019, 045. doi:10.1088/1475-7516/2019/11/045, arXiv:1910.14429.
  • Díaz Rivero et al. (2018) Díaz Rivero, A., Dvorkin, C., Cyr-Racine, F.Y., Zavala, J., Vogelsberger, M., 2018. Gravitational lensing and the power spectrum of dark matter substructure: Insights from the ETHOS N -body simulations. Phys. Rev. D 98, 103517. doi:10.1103/PhysRevD.98.103517, arXiv:1809.00004.
  • Feldmann and Spolyar (2015) Feldmann, R., Spolyar, D., 2015. Detecting dark matter substructures around the Milky Way with Gaia. MNRAS 446, 1000--1012. doi:10.1093/mnras/stu2147, arXiv:1310.2243.
  • Garrison-Kimmel et al. (2017) Garrison-Kimmel, S., Wetzel, A., Bullock, J.S., Hopkins, P.F., Boylan-Kolchin, M., Faucher-Giguère, C.A., Kereš, D., Quataert, E., Sanderson, R.E., Graus, A.S., et al., 2017. Not so lumpy after all: modelling the depletion of dark matter subhaloes by milky way-like galaxies. Monthly Notices of the Royal Astronomical Society 471, 1709--1727.
  • Gilman et al. (2019) Gilman, D., Birrer, S., Treu, T., Nierenberg, A., Benson, A., 2019. Probing dark matter structure down to 10 solar masses: flux ratio statistics in gravitational lenses with line-of-sight haloes. MNRAS 487, 5721--5738. doi:10.1093/mnras/stz1593, arXiv:1901.11031.
  • Górski et al. (2005) Górski, K.M., Hivon, E., Banday, A.J., Wandelt, B.D., Hansen, F.K., Reinecke, M., Bartelmann, M., 2005. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. The Astrophysical Journal 622, 759--771. doi:10.1086/427976, arXiv:arXiv:astro-ph/0409513.
  • Helmi et al. (2019) Helmi, A., Irwin, M., Deason, A., Balbinot, E., Belokurov, V., Bland-Hawthorn, J., Christlieb, N., Cioni, M.R.L., Feltzing, S., Grebel, E.K., Kordopatis, G., Starkenburg, E., Walton, N., Worley, C.C., 2019. 4MOST Consortium Survey 1: The Milky Way Halo Low-Resolution Survey. The Messenger 175, 23--25. doi:10.18727/0722-6691/5120, arXiv:1903.02467.
  • Hezaveh et al. (2016) Hezaveh, Y., Dalal, N., Holder, G., Kisner, T., Kuhlen, M., Perreault Levasseur, L., 2016. Measuring the power spectrum of dark matter substructure using strong gravitational lensing. JCAP 2016, 048. doi:10.1088/1475-7516/2016/11/048, arXiv:1403.2720.
  • Hopkins (2015) Hopkins, P.F., 2015. A new class of accurate, mesh-free hydrodynamic simulation methods. MNRAS 450, 53--110. doi:10.1093/mnras/stv195, arXiv:1409.7395.
  • Hopkins et al. (2018) Hopkins, P.F., Wetzel, A., Kereš, D., Faucher-Giguère, C.A., Quataert, E., Boylan-Kolchin, M., Murray, N., Hayward, C.C., Garrison-Kimmel, S., Hummels, C., et al., 2018. Fire-2 simulations: physics versus numerics in galaxy formation. Monthly Notices of the Royal Astronomical Society 480, 800--863.
  • Ibata et al. (2002) Ibata, R.A., Lewis, G.F., Irwin, M.J., Quinn, T., 2002. Uncovering cold dark matter halo substructure with tidal streams. MNRAS 332, 915--920. doi:10.1046/j.1365-8711.2002.05358.x, arXiv:astro-ph/0110690.
  • Karukes et al. (2020) Karukes, E.V., Benito, M., Iocco, F., Trotta, R., Geringer-Sameth, A., 2020. A robust estimate of the Milky Way mass from rotation curve data. JCAP 2020, 033. doi:10.1088/1475-7516/2020/05/033, arXiv:1912.04296.
  • Kingma and Ba (2014) Kingma, D.P., Ba, J., 2014. Adam: A Method for Stochastic Optimization. arXiv e-prints , arXiv:1412.6980arXiv:1412.6980.
  • Kipper et al. (2021) Kipper, R., Tenjes, P., Tempel, E., de Propris, R., 2021. Non-equilibrium in the solar neighbourhood using dynamical modelling with Gaia DR2. MNRAS 506, 5559--5572. doi:10.1093/mnras/stab2104, arXiv:2106.07076.
  • Kipper et al. (2020) Kipper, R., Tenjes, P., Tuvikene, T., Ganeshaiah Veena, P., Tempel, E., 2020. Quantifying torque from the milky way bar using gaia dr2. Monthly Notices of the Royal Astronomical Society 494, 3358--3367.
  • Knollmann and Knebe (2009) Knollmann, S.R., Knebe, A., 2009. Ahf: Amiga’s halo finder. The Astrophysical Journal Supplement Series 182, 608.
  • Lin et al. (2017) Lin, T., Goyal, P., Girshick, R.B., He, K., Dollár, P., 2017. Focal loss for dense object detection. CoRR abs/1708.02002. URL:, arXiv:1708.02002.
  • Mikuni and Canelli (2021) Mikuni, V., Canelli, F., 2021. Unsupervised clustering for collider physics. Physical Review D 103, 092007.
  • Mirabal and Bonaca (2021) Mirabal, N., Bonaca, A., 2021. Machine-learned dark matter subhalo candidates in the 4FGL-DR2: search for the perturber of the GD-1 stream. JCAP 2021, 033. doi:10.1088/1475-7516/2021/11/033, arXiv:2105.12131.
  • Moliné et al. (2017) Moliné, Á., Sánchez-Conde, M.A., Palomares-Ruiz, S., Prada, F., 2017. Characterization of subhalo structural properties and implications for dark matter annihilation signals. MNRAS 466, 4974--4990. doi:10.1093/mnras/stx026, arXiv:1603.04057.
  • Necib et al. (2020) Necib, L., Ostdiek, B., Lisanti, M., Cohen, T., Freytsis, M., Garrison-Kimmel, S., Hopkins, P.F., Wetzel, A., Sanderson, R., 2020. Evidence for a vast prograde stellar stream in the solar vicinity. Nature Astronomy 4, 1078--1083. doi:10.1038/s41550-020-1131-2, arXiv:1907.07190.
  • Neyman and Pearson (1992) Neyman, J., Pearson, E.S., 1992. On the problem of the most efficient tests of statistical hypotheses, in: Breakthroughs in statistics. Springer, pp. 73--108.
  • Oñorbe et al. (2014) Oñorbe, J., Garrison-Kimmel, S., Maller, A.H., Bullock, J.S., Rocha, M., Hahn, O., 2014. How to zoom: bias, contamination and Lagrange volumes in multimass cosmological simulations. MNRAS 437, 1894--1908. doi:10.1093/mnras/stt2020, arXiv:1305.6923.
  • Ostdiek et al. (2020) Ostdiek, B., Necib, L., Cohen, T., Freytsis, M., Lisanti, M., Garrison-Kimmmel, S., Wetzel, A., Sanderson, R.E., Hopkins, P.F., 2020. Cataloging accreted stars within Gaia DR2 using deep learning. Astron. Astrophys. 636, A75. doi:10.1051/0004-6361/201936866, arXiv:1907.06652.
  • Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., Ashdown, M., Aumont, J., Baccigalupi, C., Ballardini, M., Banday, A.J., Barreiro, R.B., Bartolo, N., Basak, S., Battye, R., Benabed, K., Bernard, J.P., Bersanelli, M., Bielewicz, P., Bock, J.J., Bond, J.R., Borrill, J., Bouchet, F.R., Boulanger, F., Bucher, M., Burigana, C., Butler, R.C., Calabrese, E., Cardoso, J.F., Carron, J., Challinor, A., Chiang, H.C., Chluba, J., Colombo, L.P.L., Combet, C., Contreras, D., Crill, B.P., Cuttaia, F., de Bernardis, P., de Zotti, G., Delabrouille, J., Delouis, J.M., Di Valentino, E., Diego, J.M., Doré, O., Douspis, M., Ducout, A., Dupac, X., Dusini, S., Efstathiou, G., Elsner, F., Enßlin, T.A., Eriksen, H.K., Fantaye, Y., Farhang, M., Fergusson, J., Fernandez-Cobos, R., Finelli, F., Forastieri, F., Frailis, M., Fraisse, A.A., Franceschi, E., Frolov, A., Galeotta, S., Galli, S., Ganga, K., Génova-Santos, R.T., Gerbino, M., Ghosh, T., González-Nuevo, J., Górski, K.M., Gratton, S., Gruppuso, A., Gudmundsson, J.E., Hamann, J., Handley, W., Hansen, F.K., Herranz, D., Hildebrandt, S.R., Hivon, E., Huang, Z., Jaffe, A.H., Jones, W.C., Karakci, A., Keihänen, E., Keskitalo, R., Kiiveri, K., Kim, J., Kisner, T.S., Knox, L., Krachmalnicoff, N., Kunz, M., Kurki-Suonio, H., Lagache, G., Lamarre, J.M., Lasenby, A., Lattanzi, M., Lawrence, C.R., Le Jeune, M., Lemos, P., Lesgourgues, J., Levrier, F., Lewis, A., Liguori, M., Lilje, P.B., Lilley, M., Lindholm, V., López-Caniego, M., Lubin, P.M., Ma, Y.Z., Macías-Pérez, J.F., Maggio, G., Maino, D., Mandolesi, N., Mangilli, A., Marcos-Caballero, A., Maris, M., Martin, P.G., Martinelli, M., Martínez-González, E., Matarrese, S., Mauri, N., McEwen, J.D., Meinhold, P.R., Melchiorri, A., Mennella, A., Migliaccio, M., Millea, M., Mitra, S., Miville-Deschênes, M.A., Molinari, D., Montier, L., Morgante, G., Moss, A., Natoli, P., Nørgaard-Nielsen, H.U., Pagano, L., Paoletti, D., Partridge, B., Patanchon, G., Peiris, H.V., Perrotta, F., Pettorino, V., Piacentini, F., Polastri, L., Polenta, G., Puget, J.L., Rachen, J.P., Reinecke, M., Remazeilles, M., Renzi, A., Rocha, G., Rosset, C., Roudier, G., Rubiño-Martín, J.A., Ruiz-Granados, B., Salvati, L., Sandri, M., Savelainen, M., Scott, D., Shellard, E.P.S., Sirignano, C., Sirri, G., Spencer, L.D., Sunyaev, R., Suur-Uski, A.S., Tauber, J.A., Tavagnacco, D., Tenti, M., Toffolatti, L., Tomasi, M., Trombetti, T., Valenziano, L., Valiviita, J., Van Tent, B., Vibert, L., Vielva, P., Villa, F., Vittorio, N., Wandelt, B.D., Wehus, I.K., White, M., White, S.D.M., Zacchei, A., Zonca, A., 2020. Planck 2018 results. VI. Cosmological parameters. Astron. Astrophys. 641, A6. doi:10.1051/0004-6361/201833910, arXiv:1807.06209.
  • Sakurada and Yairi (2014) Sakurada, M., Yairi, T., 2014. Anomaly detection using autoencoders with nonlinear dimensionality reduction, in: MLSDA’14.
  • Sanderson et al. (2020) Sanderson, R.E., Wetzel, A., Loebman, S., Sharma, S., Hopkins, P.F., Garrison-Kimmel, S., Faucher-Giguère, C.A., Kereš, D., Quataert, E., 2020. Synthetic Gaia Surveys from the FIRE Cosmological Simulations of Milky Way-mass Galaxies. The Astrophysical Journal Supp. Series 246, 6. doi:10.3847/1538-4365/ab5b9d, arXiv:1806.10564.
  • Schneider et al. (2012) Schneider, A., Smith, R.E., Macciò, A.V., Moore, B., 2012. Non-linear evolution of cosmological structures in warm dark matter models. MNRAS 424, 684--698. doi:10.1111/j.1365-2966.2012.21252.x, arXiv:1112.0330.
  • Shen et al. (2021) Shen, J., Eadie, G.M., Murray, N., Zaritsky, D., Speagle, J.S., Ting, Y.S., Conroy, C., Cargile, P.A., Johnson, B.D., Naidu, R.P., Han, J.J., 2021. The Mass of the Milky Way from the H3 Survey. arXiv e-prints , arXiv:2111.09327arXiv:2111.09327.
  • Shih et al. (2021) Shih, D., Buckley, M.R., Necib, L., Tamanas, J., 2021. Via Machinae: Searching for Stellar Streams using Unsupervised Machine Learning. MNRAS doi:10.1093/mnras/stab3372, arXiv:2104.12789.
  • Shih et al. (2021) Shih, D., Buckley, M.R., Necib, L., Tamanas, J., 2021. Via Machinae: Searching for Stellar Streams using Unsupervised Machine Learning. Monthly Notices of the Royal Astronomical Society URL:, doi:10.1093/mnras/stab3372, arXiv: stab3372.
  • Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R., 2014. Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research 15, 1929--1958.
  • Van Tilburg et al. (2018) Van Tilburg, K., Taki, A.M., Weiner, N., 2018. Halometry from astrometry. JCAP 2018, 041. doi:10.1088/1475-7516/2018/07/041, arXiv:1804.01991.
  • Vattis et al. (2020) Vattis, K., Toomey, M.W., Koushiappas, S.M., 2020. Deep learning the astrometric signature of dark matter substructure. arXiv e-prints , arXiv:2008.11577arXiv:2008.11577.
  • Vogelsberger et al. (2016) Vogelsberger, M., Zavala, J., Cyr-Racine, F.Y., Pfrommer, C., Bringmann, T., Sigurdson, K., 2016. ETHOS - an effective theory of structure formation: dark matter physics as a possible explanation of the small-scale CDM problems. MNRAS 460, 1399--1416. doi:10.1093/mnras/stw1076, arXiv:1512.05349.
  • Wang et al. (2020a) Wang, J., Bose, S., Frenk, C.S., Gao, L., Jenkins, A., Springel, V., White, S.D.M., 2020a. Universal structure of dark matter haloes over a mass range of 20 orders of magnitude. Nature 585, 39--42. doi:10.1038/s41586-020-2642-9, arXiv:1911.09720.
  • Wang et al. (2020b) Wang, W., Han, J., Cautun, M., Li, Z., Ishigaki, M.N., 2020b. The mass of our Milky Way. Science China Physics, Mechanics, and Astronomy 63, 109801. doi:10.1007/s11433-019-1541-6, arXiv:1912.02599.
  • Wetzel et al. (2016) Wetzel, A.R., Hopkins, P.F., Kim, J.h., Faucher-Giguère, C.A., Kereš, D., Quataert, E., 2016. Reconciling Dwarf Galaxies with CDM Cosmology: Simulating a Realistic Population of Satellites around a Milky Way-mass Galaxy. The Astrophysical Journal Letters 827, L23. doi:10.3847/2041-8205/827/2/L23, arXiv:1602.05957.
  • Yoon et al. (2011) Yoon, J.H., Johnston, K.V., Hogg, D.W., 2011. Clumpy Streams from Clumpy Halos: Detecting Missing Satellites with Cold Stellar Structures. The Astrophysical Journal 731, 58. doi:10.1088/0004-637X/731/1/58, arXiv:1012.2884.
  • Zavala and Frenk (2019) Zavala, J., Frenk, C.S., 2019. Dark Matter Haloes and Subhaloes. Galaxies 7, 81. doi:10.3390/galaxies7040081, arXiv:1907.11775.
  • Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., Reinecke, M., Rosset, C., Hivon, E., Gorski, K., 2019. healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in python.

    Journal of Open Source Software 4, 1298.

    URL:, doi:10.21105/joss.01298.
  • Zybin et al. (1999) Zybin, K.P., Vysotsky, M.I., Gurevich, A.V., 1999. The fluctuation spectrum cut-off in a neutralino dark matter scenario. Physics Letters A 260, 262--268. doi:10.1016/S0375-9601(99)00434-X.