The Spur and the Gap in GD-1: Dynamical Evidence for a Dark Substructure in the Milky Way Halo

Unlike the observed GD-1, a globular cluster disrupting on the GD-1 orbit in a simple—analytic and smooth—galaxy creates a stream that is also smooth (Price-Whelan & Bonaca 2018). This model follows stars as they leave the progenitor and accounts for their epicylic motion relative to the progenitor's orbit (Küpper et al. 2008, 2010; Fardal et al. 2015). The resulting pattern of over- and underdensities is much more uniform than the observed stream, so the full extent of density variations in GD-1 cannot be simply explained by the process of globular cluster disruption alone. As inhomogeneities can also be introduced into a stellar stream by adding a perturbation to the gravitational potential (e.g., Siegal-Gaskins & Valluri 2008), in this section, we present a model of the GD-1 stream that had a recent close encounter with a dense, massive object.

As a first step in creating a model of the GD-1 stream, we follow Price-Whelan & Bonaca (2018) in finding the orbit of the GD-1 progenitor by fitting the six-dimensional phase-space distribution of GD-1 stars. We assume a spherical logarithmic potential with a circular velocity of 225 km s⁻¹ for orbit integration. This simple gravitational potential is very close to the best-fit model of GD-1 (Koposov et al. 2010; Bowden et al. 2015), and it also allows much faster force evaluations than the standard multicomponent model of the Milky Way. The present-day location of the GD-1 progenitor is not firmly established (e.g., Price-Whelan & Bonaca 2018; Webb & Bovy 2019), so we adopt a scenario in which the progenitor is completely dissolved because it qualitatively reproduces a number of global GD-1 features (see Price-Whelan & Bonaca 2018).

We assume that the GD-1 progenitor was a globular cluster of initial mass 7 × 10⁴ M_⊙. In our model, it started losing stars through evaporation 3 Gyr ago and completely disrupted 300 Myr before the present day. We follow the progenitor's dissolution by releasing test particles from its Lagrange points (the mean ejection radius is ∼150 pc initially) and produce a model of the stream following Fardal et al. (2015). In this approach, the tidal radius from which the stars are being ejected is a function of the progenitor's mass and position in the Galaxy. Although idealized, such models capture detailed properties of the more realistic N-body simulations of disrupting globular clusters (Küpper et al. 2012). The present-day distribution of test particles is shown in GD-1 coordinates in the bottom panel of Figure 1. Had the progenitor survived to the present, it would be located at ϕ₁ = −20°. Instead, this model has a gap at that location, which coincides with the G-20 gap observed in GD-1. The progenitor's initial mass and time of disruption were chosen to reproduce the stream width and the morphology of the more depleted observed gap.

However, the observed GD-1 stream has two prominent gaps. To produce a model stream that also has a gap coinciding with the observed G-40 gap, our model also includes a massive and dense object on an orbit that crosses GD-1. The parameters of the encounter were chosen to reproduce the observed morphology of the gap-and-spur feature in GD-1. The perturber is represented by a Hernquist (1990) sphere of mass 5 × 10⁶ M_⊙ and scale radius 10 pc. Its closest approach to GD-1 happened 495 Myr ago with a relative distance of 15 pc, which is smaller than the stream width. The perturber came closest to GD-1 stars that are presently at ϕ₁ ≈ −38°. During the encounter, nearby stars received a velocity kick from the perturber and started moving toward the location of its closest approach. In case of a weaker perturbation, e.g., one produced by a more diffuse perturber, the most significant component of the velocity kick is along the stream, which changes the orbital period of the affected stars (Erkal & Belokurov 2015a). On one side of the perturber, the affected stars have shorter orbital periods and hence speed by the unaffected stars, while on the other side, they take longer to orbit the Galaxy and lag behind the unaffected stream stars. This creates a gap at the projected location of the closest approach, with a pileup of stars on either side of the gap creating a signature double-horned profile (Carlberg 2012). However, the perturber in our model is dense, so it imparts a significant velocity kick perpendicular to the stream, as well as along the stream. This leads to a loop of stars straying beyond the unperturbed stream track. At the present, this loop is viewed nearly edge-on and looks like the observed spur (Figure 1, bottom).

The stream model in the bottom panel of Figure 1 qualitatively matches many features in the observed GD-1 stream (Figure 1, top). Not only are both of the most prominent gaps reproduced at the right location and with the right size, but their density contrast is matched as well. The G-20 gap, modeled as a disrupted progenitor, is almost completely depleted, while the G-40 gap, the location of the impact, still retains some stars. Furthermore, the model features a spur of the correct offset from the main stream and correct length. It is not a perfect model; for example, the model stream extends past the observed extent of GD-1. Still, this model is a remarkably realistic rendition of the observed GD-1, so we next quantitatively explore the range of impact parameters that produce a good match to the observed stream.

In our fiducial model of the GD-1 stellar stream (Section 3.1), the encounter with a perturber introduced a gap in the stream and ejected a spur of GD-1 stars beyond the main stream. In this section, we constrain the range in the perturber's mass and size; its impact parameter, velocity, and time of perturbation that reproduce well the location and width of the gap; and the location and extent of the spur.

To efficiently explore the allowed parameter space, in this section, we assume that GD-1 stream stars share the orbit with their progenitor. This allows for even faster generation of a stream model, at the expense of less realistic density contrasts along the stream. In the third row of Figure 2, we show our fiducial model of perturbed GD-1 under this assumption. The perturber parameters have the same normalization as the streak-line model presented in the previous section, but the encounter velocity and impact parameter angles are different. This modification was necessary because in general, stream tracks are offset from their progenitors' orbits (Sanders & Binney 2013). Despite this difference, modeling the GD-1 stream with a single orbit results in features similar to those seen in the streak-line model. This simpler model also features a loop of stars removed from their original orbit, which projects to the GD-1 spur location at ϕ₁ ≳ −40°, opens a gap at ϕ₁ ≈ −40°, and reconnects to the leading tail at ϕ₁ ≲ −40°.

Zoom In Zoom Out Reset image size

With a method at hand to quickly generate stream models that reproduce the basic features observed in GD-1, we explore how the stream morphology depends on the impactor's properties. We consider five parameters: the time of impact, T, and the perturber's mass (M) and its impact parameter (b), velocity (V), and size (r_s). Each column of Figure 2 shows models where one of the parameters is changed, while the others are kept at their fiducial values. Parameter values increase from top to bottom (as labeled at the top left of every panel), with the fiducial values shown in the third row. Most of the presented models preserve the angle between the loop and the unperturbed stream, while the main differences are in the width of the gap and the extent of the loop.

As the mass of the perturber increases, the imparted velocity kicks to the stream stars are larger, and the resulting loop and gap also increase in size. The effects of other parameters are similar; for example, increasing the perturber's size (making it less dense) decreases the loop and gap sizes. Indeed, the bottom right panel of Figure 2 shows a model perturbed with an object following the ΛCDM concentration–mass relation (Diemer & Joyce 2019), and it shows no signature of a loop. Given time, both the loop and the gap grow in size. However, older loops are more aligned with the stream (Figure 2, bottom left panel) and hence more difficult to detect observationally. The presence of an observable spur alone already implies that GD-1 had a recent encounter with a dense perturber.

Decreasing the velocity of the perturber (which also increases the velocity kicks) produces larger loops and gaps, and this effect is nearly degenerate with an increase in the perturber's mass at sufficiently large velocities. However, at very low velocities, the whole loop morphology is changed (top panel in fourth column, Figure 2), likely because the encounter is no longer impulsive. This means that the observed spur morphology cannot be explained by an object of an arbitrarily low mass moving slowly. Likewise, decreasing the impact parameter has a unique signature too: the resulting loop is larger, but the gap size remains constant. Different dependences of the loop and gap properties on the perturber's mass and its impact parameter suggest that both parameters can be constrained in the case of GD-1, where we observe both features.

Different encounter parameters produce a unique impact on the GD-1 morphology, so we next search this space for parameter combinations allowed by the data. Ideally, we would like to create a model of the stream and compare it to the data directly. However, the adopted method for generating stream models is not sufficiently realistic for full forward modeling. Instead, we devised a set of criteria that allow us to compare whether a conceptual stream model represents the observed structure of GD-1.

First, we compare ϕ₂ positions perpendicular to the stream of perturbed model stars to the position of observed stars in the spur, defined as a spline that goes through control points of the stream between ϕ₁ = −50° and −39° and control points on the spur between ϕ₁ = −39° and −30°. The positions of likely GD-1 members and this spline are shown in the top left panel of Figure 3 as black points and a gray line, respectively. The positions of stars in the fiducial model (A) and two other models (B and C) are shown in the bottom rows, with darker colors marking stars with a significant change in energy due to the encounter (bracketed by the third percentile on the trailing end and the 92nd percentile of the energy difference on the leading end of the stream). Quantitatively, we define the spur χ²,

$$\begin{eqnarray}&&{\chi }_{\mathrm{spur}}^{2}=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{{\phi }_{2,i}-{\phi }_{2,\mathrm{spline}}({\phi }_{1,i})}{{\sigma }_{\mathrm{spur}}}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

for N model stars with positions (ϕ_1,i, ϕ_2,i) and the adopted width of the spur as σ_spur = 02. To ensure that the model spur is long enough, we only consider models where at least one star has ϕ₁ > −32° and ϕ₂ > 08.

Zoom In Zoom Out Reset image size

Next, we compare the location and width of the observed and model gaps. The observed gap profile is shown in the top panel of the middle column of Figure 3 (black points) and is well represented by a top-hat profile centered on ϕ_1,gap = −405 and w_gap = 85 wide (gray line). The gap profiles of models A–C have similar positions and widths; however, the density contrast between the gap and the spur is much larger (bottom rows). As the predictions of our stream models regarding density are simplistic, we define the gap χ² as

$$\begin{eqnarray}&&{\chi }_{\mathrm{gap}}^{2}=\displaystyle \frac{1}{{N}_{\mathrm{bin}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{bin}}}{\left(\displaystyle \frac{{N}_{\mathrm{model}}({\phi }_{1,i})-{N}_{\mathrm{top}}({\phi }_{1,i})}{{\sigma }_{\mathrm{gap},i}}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where i denotes N_bin = 29 bins in the range −60° < ϕ₁ < −20°. Here N_model(ϕ_1,i) is the number of model stars in bin i, and σ_gap,i is the associated Poisson uncertainty. Here ${N}_{\mathrm{top}}({\phi }_{1,i})\equiv {N}_{\mathrm{top}}({n}_{\mathrm{base}},{n}_{\mathrm{hat}},{\phi }_{1,\mathrm{gap}},{w}_{\mathrm{gap}}| {\phi }_{1,i})$ is the number of stars expected in bin i from a top-hat distribution with the position ϕ_1,gap and width w_gap adopted from the observed profile but with the base level, n_base, given by the median of the model profile outside of the gap (−55° < ϕ₁ < −45° and −35° < ϕ₁ < −25°), and n_hat is the median of the model profile in the gap (−43° < ϕ₁ < −37°). We further require the density contrast n_hat/n_base to be at least 0.5.

We combine the spur-and-gap contributions to the surrogate log likelihood as $\mathrm{ln}{ \mathcal L }=-({\chi }_{\mathrm{spur}}^{2}+{\chi }_{\mathrm{gap}}^{2})$. Even though this likelihood is an approximation to the formal likelihood, which would compare the positions of model stream stars to the observed ones, it is expected to favor models that reproduce well the features seen in the data (spur and gap). Therefore, we use an ensemble MCMC sampler (Foreman-Mackey et al. 2013) to find the allowed range in the parameters of interest—the perturber's mass, size, velocity, impact parameter, and impact time—while marginalizing over the orientation angles and impact location. We started 200 walkers and advanced them for 5000 steps, keeping the last 2000 steps for analysis. In principle, the resulting chain provides posterior samples, but since this is a highly idealized search of the parameter space, we only provide plausible ranges of parameters, instead of showing their two-dimensional distributions. Specifically, in each panel of the corner plot (Figure 3, right), we show the two-dimensional convex hull of all models with a likelihood above the fifth percentile (gray shaded regions). Therefore, Figure 3 contains likelihood, rather than posterior, information. Applied to the fiducial model of GD-1 (presented in Figure 1), this method recovers the true encounter parameters.

High-likelihood models of the GD-1 stream have encounter parameters expected from the stream's sensitivity to different parameters (explored in Figure 2). Recent encounters are favored, with most of the models having been perturbed within the last 1 Gyr and none more than 2 Gyr ago. The perturbing object itself is massive ($5.5\lesssim \mathrm{log}M/{M}_{\odot }\lesssim 8$) and dense (r_s ≲ 20 pc). A range of velocities are allowed, as the data dismiss only the slowest-moving perturbers (V ≳ 50 km s⁻¹). The closest approach was extremely close to the stream, with the impact parameter smaller than b ≲ 50 pc.

Many of the inferred encounter parameters are correlated. For example, more massive perturbers allow for larger impact parameters and scale radii and older encounters. In addition, the current data place very weak constraints on the perturber's velocity. This is hardly surprising given the similar effect these parameters have on the appearance of the spur and gap (Figure 2). In the next section, we discuss new observables that can further constrain the GD-1 perturber.

The perturber properties presented in the previous section were constrained by the spatial distribution of the GD-1 debris alone. We now explore the kinematic signatures of the GD-1 encounter with such a perturber and discuss how to further constrain its properties with future kinematic observations.

In Figure 3, we introduce three models whose loops reproduce the observed spur-and-gap morphologies (left); consequently, their encounter parameters are in the allowed region of the parameter space (right). These particular models were selected to display the diversity in kinematic signatures of the loop allowed by the spatial data. To account for the velocity gradients present along the stream, in Figure 4, we show differences in velocity components between the perturbed and unperturbed stream for these three models. Each column displays velocity signatures of a single model, with radial velocity differences in the top row and the two proper-motion components in GD-1 coordinates in the next two rows. In all of the models, there are kinematic offsets between the loop and the unperturbed stream.

Zoom In Zoom Out Reset image size

Radial velocity signatures exhibit the largest diversity, most prominently manifested as the velocity difference between the spur and the unperturbed stream at ϕ₁ ≈ −30°. Here model A predicts that the spur moves faster than the stream, and in model B, the stream moves faster than the spur, whereas there is no difference between the two in model C. The other side of the loop has the opposite behavior, so models A and B predict bimodal radial velocity measurements in the GD-1 stream at −50° ≲ ϕ₁ ≲ −45°. Models A, B, and C have rather similar parameters overall (Figure 3, right), but the perturber's orbital plane is misaligned with that of GD-1 by 170°, 15°, and 50°, respectively. Therefore, future follow-up observations will not only further constrain all impact parameters but also provide the first strong constraints on the encounter geometry.

Velocity offsets in proper motions have the same direction across all models (Figure 4, bottom two rows). Focusing first on the ϕ₁ direction, this behavior stems from the fact that there is a gap in the stream and the stream moves in the negative ϕ₁ direction (i.e., μ_ϕ1 < 0). As Erkal & Belokurov (2015a) showed, a velocity kick along the stream (the ϕ₁ direction) makes the affected stars on the leading side of the gap speed past the unperturbed stream. And indeed, loop stars at ϕ₁ ≲ −45° have proper motion ${\mu }_{{\phi }_{1}}$ more negative than the unperturbed stream at the same ϕ₁. The opposite behavior is expected, and seen, on the other side of the loop (ϕ₁ ≳ −35°), where the affected stars lag behind the stream at less negative ${\mu }_{{\phi }_{1}}$. Details of the loop proper-motion offsets still depend on the perturber's parameters, but the direction of offsets along the stream is universal in the encounter scenario.

Because the spur is at positive ϕ₂, its proper motion perpendicular to the stream (the ϕ₂ direction) should be faster than that of the unperturbed stream stars. As expected, ${\rm{\Delta }}{\mu }_{{\phi }_{2}}\gt 0$ for all of the models at ϕ₁ ≳ −35°. The other side of the loop is at most slightly offset from the main stream track, so ${\rm{\Delta }}{\mu }_{{\phi }_{2}}$ is very small at ϕ₁ ≲ −45°. Finally, the magnitude of kinematic offsets perpendicular to the stream (${\rm{\Delta }}{\mu }_{{\phi }_{2}}$) is universally smaller than along the stream (${\rm{\Delta }}{\mu }_{{\phi }_{1}}$). This is also expected because the gap is much wider than the distance between the stream and the spur.

The encounter model makes falsifiable predictions for proper-motion kinematics in the affected region of the GD-1 stream: the spur, as a trailing part of the loop, should be moving slower than the adjacent stream stars in the ϕ₁ direction (along the stream). The predicted velocity offsets are small, but the proper motions of diffuse, cold streams have been measured at this precision with the Hubble Space Telescope (Sohn et al. 2016), so future kinematic data will rule on the perturbative origin of GD-1 features. Furthermore, the expected offsets in radial velocity are on the order of a few km s⁻¹, which is easily within reach of current spectroscopic facilities (e.g., Simon & Geha 2007). The baseline radial velocity in the perturbed part of the stream is ∼−100 km s⁻¹ (Koposov et al. 2010), allowing for confident selection of GD-1 members out of the field Milky Way population. Should the encounter scenario be confirmed, these new data will make strong predictions on the orbit and present-day position of the perturber.

A number of simplifying assumptions were made to facilitate the initial modeling of complex features observed in the GD-1 stream. While we expect these assumptions to still produce unbiased inference, there is additional information available in the present data that has yet to be employed. We will need better models to fully explain these observations, so in this section, we discuss areas for improvement in the modeling of the GD-1 system.

Our fiducial model of GD-1 is built from test particles released by the disrupting progenitor directly from the Lagrange points, instead of particles dynamically ejected from the globular cluster progenitor itself. Stream models generated under this assumption match the distribution of tidal debris (including the intrinsic density variations along the tidal tails) from direct N-body simulations while the progenitor persists (e.g., Küpper et al. 2012; Fardal et al. 2015). However, the method is yet to be tested when the progenitor is completely disrupted, such as in the case of GD-1 and most of the known halo streams. The last stages of tidal disruption can result in enhanced density variations, so to fully account for all structures observed in GD-1, we will need a full self-gravitating N-body model of the stream.

Of course, the reason we decided against employing self-gravitating models in this work is that they are computationally expensive, prohibitively so for any kind of parameter space exploration. This is why we further focused only on the perturbed region of GD-1 when inferring the properties of its perturber and assumed that all stars are on the same orbit. Reassuringly, fiducial encounter parameters produce qualitatively similar features both in an idealized stream (Figure 1) and in stars along a single orbit (Figure 2). As discussed in Section 3.2, this choice limited us to only matching positions of stream features and disregarding density along the stream. Since the density profile of the gap is also expected to contain information about the perturber (Erkal & Belokurov 2015b), going forward, we will need to have a proper generative model of the stream.

Furthermore, our treatment of the GD-1 perturber is also simplistic. Although we explored perturbers of different masses and sizes, they all follow the Hernquist (1990) density profile. This profile reduces to the point mass case in the limit of vanishing scale size and is similar to the profile of Navarro–Frenk–White (NFW; Navarro et al. 1997) dark matter halos at small radii, although at larger radii, it falls off more steeply as r⁻⁴, compared to r⁻³ for NFW halos. While the inference of the perturber's scale radius should be robust to the details of its outer density profile, future work should explore whether any observables can be traced to the perturber's density structure, as it might hold additional clues to its origin.

In the absence of a realistic stream model, our inference of the GD-1 perturber was based on a set of high-level criteria instead of directly calculating the likelihood of the observed spatial distribution of tidal debris for a given set of model parameters (e.g., likelihood developed in Bonaca et al. 2014). Spot-checking of the accepted models suggests that we erred on the conservative side and accepted a wide range of models, some of which would likely be ruled out with the full likelihood. For example, model B in Figure 3 has both sides of its loop appreciably offset from the stream track, instead of just the trailing side coinciding with the spur. Because of that, with the current approach, we only bound the allowed parameter space and remain agnostic to the relative likelihood of models within the bounds. This prohibits us from measuring probabilities associated with each model parameter. For example, in the following section, we qualitatively discuss the inferred constraints on the mass and size of the GD-1 perturber in the context of objects orbiting the Milky Way but are unable to quantify their relative probabilities. To find the preferred regions of the parameter space and recover the posterior density, future inferences will need to employ a more realistic model of GD-1 and a proper likelihood.

And finally, our search of the parameter space has not been exhaustive, so different islands may still be allowed (for example, at a lower perturber mass or older encounter). This is, of course, always true when sampling a parameter space (e.g., Hogg & Foreman-Mackey 2018). However, ours is sufficiently low-dimensional that future studies employing the formal likelihood should be able to perform a brute-force sweep of the entire parameter volume and identify all classes of perturber properties that explain GD-1 features.

Many massive objects are known or expected to orbit the Galaxy, and if any should have encountered a cold stream, the perturbation would remain recorded in the stream's density profile. So, the presence of perturbed features in GD-1 is not surprising, but their detailed structure is. In this section, we discuss what the observed structures suggest about their origin.

4.2.1. Nonimpact Scenarios

4.2.2. Luminous Objects as the GD-1 Perturber

While these noninteraction scenarios are still possible causes of structures seen in GD-1, we find the interaction scenario more plausible. The Milky Way is surrounded by ∼50 dwarf galaxies (McConnachie 2012) and ∼150 globular clusters (Harris 2010), most of which, like GD-1, reside in the Galactic halo. Thanks to the Gaia mission, a significant fraction of these satellites are now fully located in the 6D phase space (Gaia Collaboration et al. 2018b; Simon 2018; Baumgardt et al. 2019). To test whether any of these objects could have perturbed GD-1, we integrate their orbits in a fiducial Milky Way potential (Price-Whelan 2017) and show their relative distance from the GD-1 gap at ϕ₁ = −40° during the past billion yr in Figure 5. Lines are shaded according to the object's mass, with the classical dwarfs being the darkest, ultrafaints medium, and globular clusters lighter. All of them have kept at least 1 kpc away from GD-1. The relative distances shown in Figure 5 were calculated for fiducial present-day positions and velocities of satellites. Since the associated measurement uncertainties can be substantial, we also sampled the error distribution in satellite properties and measured in how many realizations any given satellite comes closer to GD-1 than the maximum permitted impact parameter of 57 pc (see Figure 3). Upon resampling the observational uncertainties, all satellites with known orbits are excluded as GD-1 perturbers at a 3σ level or higher. Orbital parameters are known for all but two globular clusters (GLIMPSE C02 and 2MASS-GC01; Baumgardt et al. 2019), both of which are faint objects located in the Galactic bulge (Bonatto & Bica 2008; Kurtev et al. 2008) and therefore unlikely GD-1 perturbers. Out of 49 dwarf galaxies present in the up-to-date catalog of McConnachie (2012), only 24 objects have their orbits determined. In summary, present data rule out all luminous satellites with known orbits as GD-1 perturbers, but follow-up spectroscopy is required to test the remaining 25 satellites.

Zoom In Zoom Out Reset image size

As it orbits the Galaxy, GD-1 crosses the disk at timescales comparable to the inferred time of perturbation (the lightest line in Figure 5 is the distance from the Galactic plane). While strong disk shocking can facilitate disruption of a diffuse globular cluster (Dehnen et al. 2004), GD-1 disk crossings are between 13 and 23 kpc from the Galactic center, where the disk density is too low to significantly impact the stream or produce sharp features such as the gap and spur. Still, giant molecular clouds (GMCs) that are orbiting in the disk plane can perturb cold stellar streams (Amorisco et al. 2016). To test whether GMCs are viable candidates for the GD-1 perturber, in Figure 6, we compare the inferred mass and size of the GD-1 perturber (light gray shaded region) to known objects in the Milky Way, including molecular clouds. Dwarf galaxies are shown as light squares (McConnachie 2012), globular clusters as medium triangles (Baumgardt & Hilker 2018), and outer-disk molecular clouds (beyond 10 kpc) as dark circles (Miville-Deschênes et al. 2017). This comparison is rather conceptual, as different classes of objects have different density profiles: for the GD-1 perturber, we show the mass and scale radius assuming a Hernquist profile; for dwarf galaxies and globular clusters, we show the total dynamical mass and half-light radius; and we show the total (high) mass and full size for molecular clouds. Keeping these caveats in mind, the most massive globular clusters and the most compact dwarf galaxies have masses and sizes comparable to the preferred values of the GD-1 perturber, but GMCs are in general too diffuse (at a given mass, their sizes are at least an order of magnitude larger than expected of the GD-1 perturber). To additionally test for an extremely dense, yet undiscovered, class of GMCs, we also created GD-1 models perturbed by a 10⁷ M_⊙ point mass moving on a circular disk orbit for the three most recent disk crossing times. These configurations result in a spur that is below the stream at ϕ₂ < 0, rather than above at ϕ₂ > 0, as observed in GD-1. Based on their low central density and expected orbits, we conclude that GMCs are unlikely to have perturbed GD-1. A dense GMC orbiting outside of the Milky Way disk at large galactocentric radii is still allowed by the data, and in Section 4.3 we discuss future tests to ascertain the nature of the GD-1 perturber.

Zoom In Zoom Out Reset image size

The prospect of ruling out known luminous objects as GD-1 perturbers based on their orbits strongly depends on the accurate knowledge of the underlying gravitational potential. In this work, we assumed a smooth and static model for the Milky Way (Price-Whelan 2017) because it reproduces well the global 6D phase-space distribution of GD-1 (Price-Whelan & Bonaca 2018). This implies that, to first order, the assumed gravitational potential is close to the effective Milky Way potential over the last 3 Gyr (the dynamical age of GD-1 in our fiducial model). Future studies will explore whether GD-1 can distinguish between different choices for a static potential (e.g., Bovy 2015). Furthermore, the presence of two massive satellites within the Milky Way, the Large Magellanic Cloud (LMC) and Sagittarius, means the potential is more complex in detail. Dynamical considerations imply that the LMC is very massive (∼10¹¹M_⊙; e.g., Kallivayalil et al. 2013; Gómez et al. 2015; Peñarrubia et al. 2016) and has already been invoked to explain deviations from the expected stream tracks of the Tucana III and Orphan streams (Erkal et al. 2018, 2019, respectively). Sagittarius is likely less massive (∼10⁹ − 10¹¹ M_⊙; Jiang & Binney 2000) but still a significant perturber in the inner Galaxy (e.g., Laporte et al. 2019a, 2019b). These massive satellites may have affected the orbit of GD-1, as well as those of luminous satellites and any dark matter subhalos. To reaffirm that the GD-1 perturber is not a luminous satellite, future work should recalculate their relative distances in a more realistic gravitational potential that includes the LMC and Sagittarius (for example, as a time-dependent expansion of basis functions; Garavito-Camargo et al. 2019).

So far, we have estimated the influence of known objects on GD-1 by directly modeling individual bodies. An alternative approach may be to consider a statistical model for the distribution function of GMCs, globular clusters, and dwarf galaxies in our Galaxy. For example, Banik & Bovy (2019) compared the power spectra of gaps in the Palomar 5 stream expected from whole populations of GMCs and dark matter subhalos. A similar analysis of GD-1 could, in principle, estimate the likelihood of the GD-1 perturber being a member of these groups. However, unlike the gaps in Palomar 5, the gap detected in GD-1 is associated with an off-stream spur. A statistical treatment of both features is beyond the scope of this work, but the framework for simultaneous modeling of a gap-and-spur feature that we developed here should provide a good starting point for such population studies in the future.

4.2.3. Dark Perturbers

Constraining the number and properties of low-mass dark matter halos is essential for understanding the nature of dark matter (Bullock & Boylan-Kolchin 2017). Interpreting the gap-and-spur feature in the GD-1 stream as a signature of a perturbation provides us a really promising candidate for a dark halo substructure. With the follow-up work outlined below, we now have a great opportunity to detect and characterize a dark matter subhalo.

The best way to test whether GD-1 features are an outcome of an encounter with a massive perturber is to obtain detailed kinematics in the perturbed region of the stream. The encounter scenario predicts velocity offsets between the spur and the stream, and they can be used to falsify this model. The expected offsets are small but measurable with a modest investment of spectroscopic time to get the radial component of the velocity (the required precision is ∼1 km s⁻¹) and a somewhat more significant astrometric commitment to measure the transverse motion precisely enough (∼0.1 mas yr⁻¹). The magnitude of these velocity offsets could already falsify other origin scenarios. For example, a velocity offset between the spur and the stream that is much larger than the velocity dispersion in a globular cluster would imply an input of energy to the system and thus rule out scenarios in which the spur is a consequence of substructure in the GD-1 progenitor. However, the ultimate test of the encounter scenario is a measurement of velocity offsets on both sides of the GD-1 gap. The offsets are predicted to monotonically increase along the loop, with the spur always lagging behind the stream in the component of the velocity along the stream (for details, see Section 3.3).

Should the encounter scenario be confirmed, details of the kinematic structure in the perturbed part of GD-1 will put very strong constraints on the orbit of the perturber (see Figure 4 for different possibilities). This would, in turn, allow the search for additional signatures of the perturber along its inferred orbit. Recently, Van Tilburg et al. (2018) showed how a massive and dense object can be detected in the halo via time-dependent weak lensing of background stars. As the magnitude of the effect depends on the object's density profile, we might directly measure the structural properties of the GD-1 perturber. Combined with a better theoretical understanding of dark matter subhalos and their density profiles, this measurement could also provide information about the particle nature of dark matter.

Ultimately, locating the perturber gravitationally would open up possibilities for electromagnetic follow-up. If the source is bright in X-rays, this might be a signature of an accreting black hole (e.g., Bailyn et al. 1995). Alternatively, a spatially coincident excess of gamma-rays might signal the annihilation of dark matter particles (similar to the results of searches at the locations of dwarf galaxies; e.g., Hooper & Linden 2015). Either way, an electromagnetic detection would better characterize the nature of the GD-1 perturber.

As discussed above, further study of the GD-1 gap-and-spur feature may provide the first opportunity to follow up an individual halo substructure. But most excitingly, these features demonstrate that cold stellar streams are extremely fine-tuned detectors, sensitive at a level that was only hoped for beforehand. In GD-1 alone, there are additional gaps that may be further evidence of gravitational perturbations (specifically, the G-20 gap is associated with a diffuse blob of GD-1 stars beyond the main stream). In addition to GD-1, there are over 40 known streams in the Milky Way halo (e.g., Grillmair & Carlin 2016; Shipp et al. 2018). In the era of Gaia, we now have both the incentive and the resources to study them all in detail. With the full network of streams, we could learn not only about individual halo substructures, but about the population as a whole.