Probability of forming gaps in the GD-1 stream by close encounters of globular clusters

One of the most intriguing properties of the GD-1 stellar stream is the existence of three gaps. If these gaps were formed by close encounters with dark matter subhalos, the GD-1 stream opens an exciting window through which we can see the size, mass, and velocity distributions of the dark matter subhalos in the Milky Way. However, in order to use the GD-1 stream as a probe of the dark matter substructure, we need to disprove that these gaps are not due to the perturbations from baryonic components of the Milky Way. Here we ran a large number of test particle simulations to investigate the probability that each of the known globular clusters (GCs) can form a GD-1-like gap, by using the kinematical data of the GD-1 stream and GCs from Gaia EDR3 and by fully taking into account the observational uncertainty. We found that the probability that all of the three gaps were formed by GCs is as low as $1.7\times10^{-5}$ and the expected number of gaps formed by GCs is only $0.057$ in our fiducial model. Our result highly disfavors a scenario in which GCs form the gaps. Given that other baryonic perturbers (e.g., giant molecular clouds) are even less likely to form a gap in the retrograde-moving GD-1 stream, we conclude that at least one of the gaps in the GD-1 stream was formed by dark matter subhalos if the gaps were formed by flyby perturbations.


INTRODUCTION
Recent large surveys have revealed dozens of stellar streams in the halo of the Milky Way (MW), which are one-dimensional substructures of halo stars resulting from disruption of merging stellar systems, such as dwarf galaxies or globular clusters (GCs) (Malhan & Ibata 2018). These stellar streams not only confirm the prediction of the standard ΛCDM cosmology that the galaxies form through hierarchical mergers, but also allow measurement of the large-and small-scale mass distribution in the Milky Way (Koposov et al. 2010;Malhan & Ibata 2019;Bonaca et al. 2019;Banik et al. 2021).
The GD-1 stream (Grillmair & Dionatos 2006) is one of the most studied stellar streams. This stream is characterized by its thinness (∼ 70 pc) and its length (more than several tens of degrees) (Carlberg & Grillmair 2013). The integrated light from the visible part of Email: dokeyuka@g.ecc.u-tokyo.ac.jp Email: khattori@ism.ac.jp Corresponding author: Yuka Doke, Kohei Hattori * These authors contributed equally to this work.
the GD-1 stream suggest that the progenitor system of this stream initially had a mass of (1.58±0.07)×10 4 M (de Boer et al. 2020), although this mass estimate might be a lower limit if the stream is much more extended than currently recognized. 1 These structural properties suggest that its progenitor system (which has been completely disrupted) was a GC-like system which began disruption at least ∼ 3 Gyr ago (Bowden et al. 2015;Erkal et al. 2016).
Recently, there have been discoveries of multiple under-dense regions or 'gaps' in the GD-1 stream (Carlberg & Grillmair 2013;de Boer et al. 2018de Boer et al. , 2020Price-Whelan & Bonaca 2018;Bonaca et al. 2019). Although there is some argument that these gaps might be due to the formation process of the stream (such as the nonuniform stripping rate; epicyclic overdensities, Küpper et al. 2010Küpper et al. , 2015 or the disruption of the progenitor star cluster within a host dwarf galaxy, Malhan et al. 2019Malhan et al. , 2021Qian et al. 2022), currently the most widely accepted hypothesis to explain these gaps is that these gaps were generated by perturbations from massive compact objects, including the dark matter subhaloes in the MW (Carlberg 2009(Carlberg , 2016Yoon et al. 2011;Erkal & Belokurov 2015a,b;Erkal et al. 2016). In general, when a stream experiences a close encounter with a dark matter subhalo with a mass of ∼ 10 5 -10 7 M , the impulsive force from the subhalo results in a differential velocity kick along the stream and forms a gap in the stream. Since the gap forming mechanism is well understood, the morphology of the gap and the velocity structure near the gap enable us to investigate the mass and size of the perturber, relative velocity of the perturber with respect to the stream, and the time of the encounter (Erkal & Belokurov 2015a). Also, the number of gaps in a stream can put a constraint on the abundance of the dark matter subhalos if all of the gaps are generated by the dark matter subhalos (Erkal & Belokurov 2015b).
To use these gaps as a probe of the dark matter, we need to disprove that these gaps are not generated by baryonic effects, such as the perturbation from the Galactic bar (Hattori et al. 2016;Price-Whelan et al. 2016), spiral arms of the MW (Banik et al. 2021), giant molecular clouds in the stellar disk (Amorisco et al. 2016), dwarf galaxies (Bonaca et al. 2019), and GCs (Bonaca et al. 2019;Banik et al. 2021). In disproving these effects, it is informative to note that the gaps are more difficult to generate if a perturber moves at a larger velocity relative to the stream; A higher-speed encounter results in a shorter interaction time with the stream and, therefore, a smaller effect. Given that the GD-1 stream is orbiting around the MW in a retrograde fashion, the perturbations from the bar, spiral arms, and giant molecular clouds -all of which show a prograde motion -have negligible effects in generating gaps in the GD-1 stream. 2 Since the dwarf galaxies are far away from the GD-1 stream at the current epoch and in the past, dwarf galaxies also have negligible effect in generating gaps in the GD-1 stream (Bonaca et al. 2019). In this regard, it is crucial to investigate whether GCs can explain the gaps in the GD-1 stream. Previously, Bonaca et al. (2019) tried to pursue this strategy by using the catalog of ∼ 150 GCs equipped with the astrometric data from Gaia DR2 (Gaia Collaboration et al. 2016. They integrated the orbit backward in time for 1 Gyr to conclude that the impacts from these GCs are negligible in the last 1 Gyr. In this paper, we updated their results by using a larger number of GCs (158 GCs) from Vasiliev & Baumgardt (2021) equipped with the astrometric data from Gaia EDR3 (Gaia Collaboration et al. 2021), and, more importantly, by adopting a long enough integration time (6 Gyr in our fiducial model). We note that Banik et al. (2021) investigated how the baryonic perturbers (including GCs, giant molecular clouds, and spiral arms) affect the density power spectrum along the GD-1 stream by running simulations of the GD-1 stream in the last 3-7 Gyr. Our work is complementary to Banik et al. (2021), because we focus on the probability that all of the three gaps in the GD-1 stream were formed by GCs.
This paper is organized as follows. In Section 2, we present the observed data of the GD-1 stream and the GCs. In Section 3, we describe how we set up our testparticle simulations. In Section 4, we show the results of our analysis, details of the gap-forming GCs, and the show-case stream models. In Section 5, we present some discussion of our results, including the estimated probability that all three gaps in the GD-1 stream were formed by GCs. In Section 6, we summarize our paper.
2. DATA Here we describe the data of the GD-1 stream and GCs. The coordinate system is defined in Appendix A.

Globular clusters
Using the astrometric data from Gaia EDR3 (Gaia Collaboration et al. 2021), Vasiliev & Baumgardt (2021) derived the position and velocity of 170 GCs in the Milky Way. We selected 158 GCs with the full 6D position-velocity data by omitting 12 GCs with incomplete data. The data include Right Ascension and Declination (α, δ), distance d, line-of-sight velocity v los , and proper motion (µ α * , µ δ ), and their associated uncertainties including the correlation between the two proper motion components. We note that the final sample of GCs contains some of the newly discovered GCs for which previous studies (e.g., Bonaca et al. 2019) have not checked whether they have experienced a close encounter with the GD-1 stream.

Candidate stars of the GD-1 stream
By using STREAMFINDER algorithm (Malhan & Ibata 2018;Malhan et al. 2018a,b), Malhan & Ibata (2019) compiled a catalog of 97 candidate stars of the GD-1 stream (see Table 2 of Malhan & Ibata 2019) for which line-of-sight velocity v los is taken from spectroscopic surveys such as SEGUE (Yanny et al. 2009) and LAMOST (Zhao et al. 2012). Since their original catalog is based on Gaia DR2 (Gaia Collaboration et al. 2018), we crossmatched these 97 stars with Gaia EDR3 (Gaia Collaboration et al. 2021) to obtain more accurate astrometric data. Judging from v los , some of their sample stars are probably not a member of the GD-1 stream (see Fig. 4d of Malhan & Ibata 2019). However, we did not make any manual selection to discard these outliers, because such a hard cut might bias our inference on the properties of the GD-1 stream. Instead, we modeled the position and velocity of the GD-1 candidate stars with a mixture model of stream stars and background stars (cf., Hogg et al. 2010).

Density along the GD-1 stream
By using the photometric data from Pan-STARRS DR1 and astrometric data from Gaia DR2, de Boer et al. (2020) derived the global properties of the GD-1 stream, such as the density, distance, and proper motions along the stream. They found three gaps along the GD-1 stream at φ 1 = −36, −20, and −3 deg in the GD-1 stream coordinate system (φ 1 , φ 2 ) (see Appendix A; Koposov et al. 2010). In this paper, we interpret that all of the three gaps are due to some sort of perturbations, although some authors interpret that the gap at φ 1 = −36 deg corresponds to the location of the alreadydisrupted progenitor (e.g., Banik et al. 2021). The main goal in this paper is to evaluate the probability that three gaps are formed within −40 deg ≤ φ 1 ≤ 0 deg by close encounters with known GCs. We note that we do not aim to re-create these three gaps rigorously at the observed locations.
3. SIMULATIONS We analyzed the close encounter with the GD-1 stream and GCs in three steps, as described in Fig. 1.
• Step 1: By using a static (unperturbed) Galactic potential, we generated (i) an unperturbed model of the GD-1 stream represented by test particles; and (ii) N MC = 1000 Monte Carlo orbit models for each GC reflecting the observed uncertainties. We treat each GC as a test particle moving in the Galactic potential in this step. • Step 2: By using the unperturbed models (i) and (ii) in Step 1, we select GC orbits that could have encountered the GD-1 stream with a small impact parameter and a small relative velocity. • Step 3: For each of the selected orbits in Step 2, we construct a time-dependent, perturbed potential, consisting of the static MW potential and the time-dependent potential from the GC that moves around the MW potential. Under this composite MW+GC potential model, we ran a test-particle simulation of the GD-1 stream perturbed by a GC due to a close encounter.
Throughout this paper, we used the AGAMA package (Vasiliev 2019) to run simulations. In the following, we describe details of these steps.

Model potential of the Milky Way
In this paper, we used the MW model potential in McMillan (2017). This model is static and axisymmetric, and it consists of atomic and molecular gas disks, thin and thick stellar disks, bulge, and a dark matter halo. In Step 3 of our experiments, we add a timedependent potential caused by the gravitational force from a GC to the above-mentioned static MW potential. In such a case, we assume that the MW potential is rigid and that the MW does not wobble, because the mass of a GC is much smaller than that of the MW.

3.2.
Step 1 (i): Unperturbed model of the GD-1 stream To generate an unperturbed model of the GD-1 stream, we need to assume the orbit of its progenitor system. We assume that the GD-1's progenitor would be at α = 148.91 deg at the current epoch if it had not been completely disrupted. This location is chosen following Banik et al. (2021). By assuming that the GD-1 stream approximately delineates the orbit of its progenitor, 3 we derived the remaining 5D information of the GD-1's progenitor at the current epoch. Specifically, we fit the 6D phase-space data of the candidates of the GD-1 stream members in Malhan & Ibata (2019) with an orbit by adequately taking into account that some fraction of the stars is non-member stars. Our best-fit orbit of the GD-1's progenitor is characterized by These quantities are broadly consistent with the result in Banik et al. (2021). In general, the orbit of the GD-1's progenitor is quite uncertain, not only because the current location (α prog , δ prog )| t=0 of the invisible (disrupted) progenitor is arbitrary, but also because the MW potential has evolved over the last several Gyr. However, we use our best-fit orbit throughout this paper, because our aim is to quantify the encounter rate of the GD-1 stream with GCs and not to reproduce the entire properties of the GD-1 stream.
Given this phase-space information, we integrated the orbit of the GD-1's progenitor backward in time from t = 0 (current epoch) to t = −T . We adopt T = 6 Gyr as the fiducial value, and thus all the figures in this paper assumes this fiducial value. We comment on how the choice of T affects our results in Section 5.2.
Under the model potential we adopted ( The position and velocity of the GD-1's progenitor at t = −T was used to create the initial condition of the GD-1 stream particles. As a simple prescription to mimic the generation of the GD-1 stream, we release 10 5 test particles at t = −T from the same position as the progenitor, (x prog , y prog , z prog )| t=−T , with a relative velocity with respect to the progenitor following an isotropic Gaussian distribution: Here, N (0, σ 2 ) represents a Gaussian distribution with mean 0 and dispersion σ 2 . We regard these 10 5 particles as the unperturbed GD-1 stream, and we integrate the orbit of these 10 5 particles forward in time from t = −T to t = 0 (current epoch) under the unperturbed Galactic potential. After some experiments, we chose σ v = (0.5 km s −1 )(T /6 Gyr) −1 so that the length of the GD-1 stream model at the current epoch is comparable to the observed extent of the GD-1 stream. We recorded the snapshot of these particles every 1 Myr and used this information in Step 2. Fig. 2 shows the stellar distribution of the unperturbed model at the current epoch. We see that the unperturbed model reproduces the observed phase-space distribution of the stars except for the apparent outlier stars. This unperturbed model is used as the benchmark model with which we quantify the strength of perturbation from GCs (see Section 3.6 and Appendix B). We note that randomly chosen 100 stars in the unperturbed model is used in Step 2 to find GCs that may have experienced a close encounter with the GD-1 stream.

3.3.
Step 1 (ii): Unperturbed orbits of the globular clusters with observational uncertainty The GCs in the catalog of Vasiliev & Baumgardt (2021) typically have ∼ 3 percent error in distance. This slight difference makes a noticeable difference in evaluating whether a given GC could have experienced a close encounter with the GD-1 stream in the last few Gyr, because a slight difference in the initial condition can cause a kpc scale difference in predicting the location of the GC in a few Gyr.
To account for the observational uncertainty, we sampled the current position and velocity of each GC for N MC = 1000 times from the error distribution. Namely, we randomly sampled (d, µ α * , µ δ , v los ) for each GC, by fully taking into account the correlation between the error in (µ α * , µ δ ). We neglected the small uncertainty in (α, δ). For each GC, we integrated the orbit backward in time from t = 0 to t = −T under the unperturbed Galactic potential. We recorded the position and velocity of each GC every 1 Myr (and used this information in Step 2). For convenience, we enumerated each of these Monte Carlo orbits by two integers (j, k), where j (0 ≤ j ≤ 169) denotes jth GC in the catalog of Vasiliev & Baumgardt (2021) and k (0 ≤ k ≤ 999) denotes kth Monte Carlo orbit. For example, (j, k) = (4, 557) corresponds to k = 557th Monte Carlo orbit of NGC1261, since NGC1261 is the 4th GC in the catalog.

3.4.
Step 2: Globular clusters with close encounters To find GCs that could have experienced a close encounter with the GD-1 stream, we used the results in Step 1 (i) and (ii). First, from the unperturbed model of the GD-1 stream computed in Step 1 (i), we randomly chose 100 particles. (We note that we sample 100 particles only once; and we use the identical 100 particles throughout Step 2.) We checked the relative distance and velocity of these 100 GD-1 particles with respect to each Monte Carlo orbit of a given GC (computed in Step 1 (ii)) as a function of time. For each Monte Carlo orbit (j, k) of the GC, we searched for a close encounter with the GD-1 stream, which we defined as a moment when at least one particle in the stream is located within d min < 0.5 kpc from the GC having a relative velocity smaller than v rel < 300 km s −1 .
As a result, we found that 1383 Monte Carlo orbits of 28 GCs have experienced a close encounter with the GD-1 stream in the last T = 6 Gyr. We note that, up  (2019) (red data points with error bar). On each panel, the the horizontal axis is φ1, while the vertical axis is either (φ2, , v los , µα * , µ δ ). As seen in the top right plot showing (φ1, v los ) distribution, some of the GD-1 candidate stars are outlier (non-member) stars with very different v los . Apart from these outlier stars, our unperturbed GD-1 stream model reproduces the global properties of the observed GD-1 stream. We note that this unperturbed GD-1 stream model is used as a reference throughout this paper.
to this step, it was unclear whether each of these GCs could form a gap, 5 which was investigated in Step 3.

3.5.
Step3: The GD-1 stream models perturbed by globular clusters For the 1383 Monte Carlo orbits (j, k) found in Step 2, we ran more detailed simulations. To efficiently run perturbed model of the GD-1 stream, we make some simplifying assumptions.
In each simulation specified by (j, k), we only consider one perturber, namely the jth GC with kth Monte Carlo orbit. In other words, we do not consider a situation where multiple perturbers exist. We treat the GD-1 stream particles as test particles that feel forces from the MW and jth GC (with kth Monte Carlo orbit). We do not take into account the self-gravity of the GD-1 stream. We assume that jth GC has a Plummer density profile with the total mass M GC,j adopted from the compilation by Holger Baumgardt. 6 We assume that the mass of GCs does not change as a function of time. We assume that the scale radius of the Plummer profile is 10 pc. We note that changing the scale radius of GCs in the range of 1-20 pc does not significantly affect our results. We assume that the GC only feels the force from the MW, and therefore the orbit of the GC in Step 2 is identical to that in Step 3. We also assume that the MW potential is rigid and does not move due to the motion of the GC, which is a natural assumption if we consider the perturbation from the GC only. 7 Under the above-mentioned assumptions, the gravitational potential that a GD-1 stream particle feels at location x and at time t is given by (10) Here, Φ MW (x) is the static MW potential (Section 3.1). The kth orbit of the jth GC is denoted as x GC,jk (t). The GC's potential is described by a Pummer potential We represent the GD-1 stream with 10 5 test particles, which is large enough to statistically robustly detect a gap in the GD-1 stream model. The initial conditions of the GD-1 stream particles are chosen in the same manner as in Step 1 (i). Importantly, we use the same Note-(a) Probability of forming a gap in the GD-1 stream. We note that P gap-forming = N gap-forming /1000. (b) Number of gap-forming orbits in our simulation. We note that we tried 1000 Monte Carlo orbits for each GC. (c) Bold-face number corresponds to the show-case example in Table 2 and Figs. 6-7.
random seed to generate the initial condition of the GD-1 stream model for all the simulations in Step 1 (i) and Step 3. We integrate the orbits of these 10 5 test particles from t = −T to t = 0 under the perturbed, timedependent potential (equation (10)). After running the simulation, we judge whether each stream model contains a GD-1-like gap or not by a simple method in Section 3.6. As a result, 57 Monte Carlo orbits of 16 GCs result in a GD-1-like gap, which will be discussed in Section 4.

Definition of a GD-1-like gap in our simulation
To assess if a given model in our simulation contains a GD-1-like gap, we compare the perturbed model in Step 3 with the unperturbed model in Step 1 (i). Because we use the same random seed to create these models, whenever we detect a notable difference between these models, the difference can be attributed to the effect of the GC perturbation.
To make a fair comparison, we first compute the linear density of the unperturbed and perturbed models of the GD-1 stream at the current epoch (t = 0) along φ 1coordinate, ρ unperturbed (φ 1 ) and ρ perturbed (φ 1 ), respectively, by using the histogram of φ 1 with a bin size of ∆φ 1 = 2 deg. We define that a GD-1-like gap is seen in the model at the location φ 1 when the perturbed model stream satisfies the following two conditions: Among the 158 GCs that we explored, we identified 16 GCs (57 Monte Carlo orbits) that formed a GD-1-like gap in the simulation. These 16 GCs are good candidates of GCs that might have formed the observed gaps in the GD-1 stream. These 16 GCs are listed in Table 1, along with the gap-forming probability P gap-forming = N gap-forming-orbit /1000, where N gap-forming-orbit denotes the number of gap-forming Monte Carlo orbits.
As we see from Table 1, P gap-forming is generally low. However, there are 6 GCs that can form a GD-1-like gap with P gap-forming ≥ 3 1000 . NGC5272 (M3) and IC4499 are especially interesting GCs, which have the highest and second highest probability (0.022 and 0.009, respectively) of forming a GD-1-like gap.
Most of the 57 Monte Carlo orbits listed in Table 1 have only one close encounter with the GD-1 stream, forming a single gap in the GD-1 stream. Intriguingly, a few Monte Carlo orbits of NGC5272 (M3) have two close encounters with the GD-1 stream, forming one visible gap and another mild under-density region.

Gap-forming probability and the GC mass
In Fig. 3(a), we show the relationship between the mass and the gap-forming probability for the 16 GCs listed in Table 1. The 6 GCs with P gap-forming ≥ 3 1000 are shown by red filled circle, while the 10 GCs with 0 < P gap-forming ≤ 2 1000 are shown by black open circle. Except for NGC5272 (M3) and IC4499, there is a mild trend that more massive GCs have a higher gap-forming probability. This trend is understandable in the following manner: Given that the maximum impact parameter to form a GD-1-like gap is larger for more massive GCs, more massive GCs have a larger chance of forming a gap. In contrast, low-mass GCs need to pass very close to the GD-1 stream to form a gap, and thus lower-mass GCs have a smaller chance of forming a gap.
The 10 GCs with 0 < P gap-forming ≤ 2 1000 have only one or two Monte Carlo orbits among 1000 trials that result in a GD-1-like gap. The mass of these GCs is 10 4.6 M M GC 10 5.6 M , which roughly covers around 20-80 percentiles of the mass of known GCs (see gray histogram in Fig. 3(b)). Thus, these 10 GCs have a typical mass of the MW GCs. Intriguingly, the two highest-P gap-forming GCs, NGC5272 (M3) and IC4499, are not very massive (10 5.6 M and 10 5.2 M , respectively). Therefore, their high value of P gap-forming is not because they are very massive but because their orbits are favorable to form a gap in the GD-1 stream.
In Fig. 3(b), we show the normalized cumulative distribution of the GC mass for GCs with different ranges of P gap-forming . As we can see from this figure, the median GC mass is approximately 10 5.0 M , 10 5.25 M , and 10 5.7 M , for GCs with P gap-forming = 0, 0 < P gap-forming ≤ 2 1000 , and 3 1000 ≤ P gap-forming , respectively. This result suggests that more massive GCs tend to have a larger probability of forming a gap in the GD-1 stream, supporting the mild trend seen in Fig. 3(a). Fig. 3(b) also shows that there is no GCs with M GC 10 4.5 M that form a gap in the GD-1 stream in our simulation. This result suggests that it is extremely difficult (P gap-forming < 1 1000 ) to form a gap in the GD-1 stream by these low-mass GCs. Fig. 4 shows the relationship between the final gap location φ 1,gap and the epoch of the encounter t encounter for the 6 GCs with P gap-forming ≥ 3 1000 . We see that all of the gap-forming encounters shown here happen at t < −1.5 Gyr. This result is consistent with the previous work by Bonaca et al. (2019) which claims that none of the GCs seems to have experienced a close encounter with the GD-1 stream in the last 1 Gyr. This result also suggests that tracing the past orbits of GCs up to  Table 1. The 6 red dots are the GCs with P gap-forming ≥ 3 1000 . (b) The normalized cumulative distribution for the mass of GCs with different ranges of P gap-forming . Those GCs with higher P gap-forming are typically more massive. Also, low-mass GCs (MGC 10 4.5 M ) can hardly form a gap. a sufficiently long time ago (up to t = −6 Gyr in our simulation) is important to claim whether or not GCs can form a visible gap in the GD-1 stream.

Epoch of the gap-forming encounters
Another intriguing result seen in Fig. 4 is that some GCs have their preferable epoch to interact with the GD-1 stream. This result is most prominently seen for NGC5272 (M3) and IC4499. For example, for all of the 22 orbits of NGC5272 (M3) that formed a GD-1-like gap, the GD-1 stream encounters this GC at t encounter −3.5 Gyr. For these 22 orbits, the location of the GD-1 gap in the φ 1 coordinate is distributed at −40 deg < φ 1 < 0 deg. Thus, an encounter with NGC5272 (M3) at t encounter −3.5 Gyr can explain any of the observed gaps at −40 deg < φ 1 < 0 deg. As another example, for all of the 9 orbits of IC4499 that formed a GD-1-like gap, the GD-1 stream encounters this GC at t encounter −1.7 Gyr. For these 9 orbits, the location of the GD-1 gap in the φ 1 coordinate is distributed at −20 deg < φ 1 < 0 deg, and none of them are distributed at φ 1 < −20 deg. Thus, an encounter with IC4499 at t encounter −1.7 Gyr can explain either of the  Table 1. The 6 GCs with P gap-forming ≥ 3 1000 are marked by various symbols. The other 10 GCs are shown by gray, filled circles. The size of the symbol represents the strength of the gap. We note that the encounter with NGC5272 (M3) at tencounter −3.5 Gyr and that with IC4499 at tencounter −1.7 Gyr can form gaps at −40 deg < φ1 < 0 deg and −20 deg < φ1 < 0 deg, respectively.
gaps at φ 1 = −3 deg or −20 deg but can hardly explain the gap at φ 1 = −36 deg. To understand these preferable epochs to form a gap, we check the orbital phase of the GD-1 stream at t −3.5 Gyr when it encounters NGC5272 (M3) and at t −1.7 Gyr when it encounters IC4499. Intriguingly, all of these encounters take place when the GD-1 stream is close to its pericenter (R 14 kpc). The GD-1 stream becomes longest near the pericentric passage, and therefore the probability of a close encounter is increased due to the enlarged 'cross section' of the GD-1 stream. As mentioned in Section 4.1, these two GCs (NGC5272 (M3) and IC4499) have the highest probability of forming a gap in our fiducial simulations. Our finding hints that, in general, GCs that can encounter the GD-1 stream near the GD-1's pericenter are a promising perturber to form a gap. Fig. 5 shows the orbital properties of the GD-1 stream and GCs investigated in this study. In both panels, the GD-1 stream is marked by blue cross and the 6 GCs with P gap-forming ≥ 3 1000 are marked with red dots. As we see in Fig. 5(a), all of the GCs with P gap-forming ≥ 3 1000 have retrograde or mildly prograde orbits (L z > −1000 kpc km s −1 ). This tendency is understandable because the GD-1 stream has a highly retrograde orbit with L z,GD-1 2800 kpc km s −1 . When a GC with a certain azimuthal angular momentum L z,GC has a close encounter with the GD-1 stream at a Galactocentric cylindrical radius R, their relative velocity is at least v rel ≥ |v φ,GD-1 − v φ,GC | = |L z,GD-1 − L z,GC |/R. The apocentric and pericentric radii. The blue cross (+) corresponds to the GD-1 stream. The red dots correspond to the 6 GCs with P gap-forming ≥ 3 1000 ; among which two GCs (NGC5272 (M3) and IC4499) with the highest-P gap-forming are highlighted with red . The black open circles correspond to the 10 GCs with 0 < P gap-forming ≤ 2 1000 . The other GCs (namely, those GCs with P gap-forming = 0) are marked with gray dots.

Gap-forming probability and the orbit of GCs
Given that the radial excursion of the GD-1 stream is 14 kpc R 20 kpc (see Fig. 5(b)), those GCs with highly prograde orbits L z,GC < −2000 kpc km s −1 have at least v rel 200 km s −1 . Thus, GCs with highly prograde orbits are hard form a gap in the GD-1 stream unless they are very massive.
We note that the GD-1 stream and NGC3201 share similar orbital properties. Indeed, Malhan et al. (2022) recently claimed that these systems are part of the same merging event dubbed 'The Arjuna/Sequoia/I'itoi merger' (see also Bonaca et al. 2021). Although their orbital similarity is intriguing, the fact that most of the high-P gap-forming GCs have very different orbital properties means that even if a GC has an orbital property similar to that of the GD-1 stream, such a GC is not necessarily a good candidate for forming a gap. Rather, as we mentioned in the previous paragraph, the relative distance and velocity at their closest approach are more important factors to form a gap.

Details on the show-case models
For an illustration purpose, for each of the 6 GCs with P gap-forming ≥ 3 1000 , we selected one show-case model. The details of the selected show-case models are summarized in Table 2. The models are named NGC2808 698, NGC3201 202, NGC5272 M3 261, IC4499 311, FSR1758 731, and NGC7089 M2 138. Here, the last three digits of these names correspond to the value of k. Figs. 6 and 7 show the current-day properties of the GD-1 stream models corresponding to these show-case models.
In Figs. 6 and 7, each show-case model is displayed with three rows. The top panel of each show-case model shows the one-dimensional density ρ(φ 1 ) with an arbitrary unit. The blue solid-line histogram shows the histogram of stars in each model. The gray dashed line shows the estimated linear density (arbitrarily scaled by a constant factor) derived in de Boer et al. (2020). The contrast between the gap and its surrounding over-dense regions in our models is similar to that in the observed GD-1 stream. Also, the widths of the gaps in our models are comparable to the observed ones.
The middle and bottom panels in each model in Figs. 6 and 7 shows the morphology of the model in (φ 1 , X) space, where X corresponds to various observables. The show-case model NGC2808 698 (Fig. 6) shows a hole-like structure in (φ 1 , φ 2 ) space, caused by the perturber (in this case NGC2808) that penetrated the stream. Interestingly, due to the hole-like structure, we see two parallel sequences of the stream at −40 deg < φ 1 < −20 deg, which is reminiscent of the observed spur-like feature in the GD-1 stream (Price-Whelan & Bonaca 2018; Bonaca et al. 2019). A similar hole-like structure in (φ 1 , φ 2 ) space is also seen in the show-case model NGC3201 202. In this case, a hole-like feature is also seen in (φ 1 , d) space, indicating that this feature is a three-dimensional structure. Given that we see a hole-like structure in multiple models, it may be one of the generic features that GCs can form.
Among the 6 GCs, NGC2808, FSR1758, and NGC7089 (M2) are the most massive GCs with M > 6 × 10 5 M . Figs. 6 and 7 show that these massive GCs can form a prominent gap, even if the relative velocity of the encounter is as large as 300 km s −1 (see also Table. 2). In contrast, the show-case model IC4499 311 results in a clear but narrow gap, due to (i) the relatively small mass of IC4499 (1.55 × 10 5 M ); and (ii) the relatively recent encounter (∼ 1.7 Gyr ago). The 16 GCs with P gap-forming ≥ 1 1000 can possibly form a gap in the GD-1 stream. If we assume, as a working hypothesis, that all three gaps in the GD-1 stream were created by the perturbation from these GCs, we can estimate its probability by using P gap-forming listed in Table 1: P (GCs formed 3 gaps in the GD-1 stream) = l∈{4,9,...,165} m > l n > m [P gap-forming (l) Here, P gap-forming (j) corresponds to the gap-forming probability of jth GC, and {4, 9, . . . , 165} corresponds to the set of j listed in Table 1. Given this tiny probability, our results suggest that perturbations from GCs are difficult to explain all three gaps in the GD-1 stream. Because other baryonic effects (e.g., from spiral arms, the Galactic bar, giant molecular clouds, or dwarf galaxies) are even more unlikely to form a gap (see Section 1), our results favor a scenario in which at least one of the gaps in the GD-1 stream were formed by dark matter subhalos (Carlberg 2009(Carlberg , 2016Erkal & Belokurov 2015a,b;Erkal et al. 2016;Bonaca et al. 2019;Banik et al. 2021).
Our results can also be used to estimate the expected number of gaps formed by GCs in the GD-1 stream, N gap (GCs). By assuming that a GC can form at most one gap (and no GCs can form multiple gaps), we have N gap (GCs) = l∈{4,9,...,165} P gap-forming (l) = 0.057. (13) We note that the gaps considered in this paper satisfy condition (A) in Section 3.6, namely ρ perturbed (φ1) ρ unperturbed (φ1) < 0.8. As a reference, Erkal et al. (2016) estimated the number of gaps in the GD-1 stream formed by dark matter subhalos, N gap (subhalos). According to their Table 2, the expected number of gaps with ρ perturbed (φ1) ρ unperturbed (φ1) < 0.75 formed by dark matter subhalos with (10 5 -10 9 )M is N gap (subhalos) = 0.6, which is ∼ 10 times larger than our estimate of N gap (GCs). Although their value of N gap (subhalos) is still smaller than 3 (the observed number of gaps), this comparison also favors dark matter subhalos as the cause of the GD-1's gaps.

Dynamical age of the GD-1 stream
In this paper, we assume that all the stars in the GD-1 stream escaped from the progenitor system at t = −T . The dynamical age of the stream, T , is set to be T = 6 Gyr in the main analysis of this paper. To check how our choice of T affects our result, we ran additional simulations with T = 1, 2, 3, and 4 Gyr (see Appendix C). (We note that we kept our description in Section 3 as general as possible so that the readers can see how the change in T affects the numerical setup of the simulations.) For simulations with T = 1 Gyr and T = 2 Gyr,   we found that all the GCs have P gap-forming = 0. For simulations with T = 3 Gyr, we found only two GCs have P gap-forming ≥ 1 1000 (see Table. 3), which means that GCs can form at most two gaps. For simulations with T = 4 Gyr, we found 5 GCs with P gap-forming ≥ 1 1000 (see Table. 3). By combining the result from our fiducial models with T = 6 Gyr, we obtain P (GCs formed 3 gaps in GD-1) and We see that both the total probability and the expected number of gaps become smaller if we adopt a smaller value of T . We can understand this tendency in two ways. First, younger streams have fewer opportunities to interact with GCs. Second, it takes some time for a gap to grow and become visible. We note that the previous work by Bonaca et al. (2019) found that no GCs experienced a close encounter with the GD-1 stream in the last 1 Gyr. Their result is consistent with our result with T = 1 Gyr. Our results suggest that adopting a longer integration time increases the chance that GCs can form a gap in the GD-1 stream, but it is extremely hard to explain three gaps only by the GCs, even if we adopt a long integration time of T = 6 Gyr. We note that Banik et al. (2021) investigated the past orbit of the GD-1 stream up to t = −7 Gyr. However, they aimed to assess the power spectrum of the density along the GD-1 steam; and not to focus on the individual gaps. Thus, our work is complementary to their work.

Choice of the model Galactic potential
In this paper, we used a model MW potential in McMillan (2017) as the fiducial model. In order to check how our results are affected by the chosen MW potential, we did the same simulations with T = 6 Gyr but with potential models in Bovy (2015) and Piffl et al. (2014). As a result, we found no dramatic changes from our fiducial simulations in the gap-forming probability P (GCs formed 3 gaps in GD-1) We note that the list of GCs that form the gaps (i.e., the list of GCs in Table 1 in our fiducial model) slightly changes if we adopt different MW potentials. However, some GCs seem to be more likely to form gaps. For example, if we adopt the Galactic potential model in Bovy (2015), the three important GCs with the highest P gap-forming are IC4499, NGC7089 (M2), and NGC3201, all of which appear in Table 1. Also, if we adopt the Galactic potential model in Piffl et al. (2014), the two important GCs with the highest P gap-forming are IC4499 and NGC6101, both of which appear in Table 1. Intriguingly, the gap-forming probability of IC4499 is always high P gap-forming 0.01, independent of the adopted potential.

Caveats in our analysis
As discussed in Section 5.1, we found that it is extremely unlikely that all three gaps in the GD-1 stream are formed by known GCs. Here we discuss the limitation of our analysis and possible future directions.
In our simulation, we treated the stars in the GD-1 stream as test particles that feel the gravitational force from the MW and a perturbing GC. We assume that all the stars were stripped from the center of the GClike progenitor system at t = −T . In reality, the stripping process may be continuous, and the stars escape from the progenitor from the inner and outer Lagrange points (Küpper et al. 2010(Küpper et al. , 2012(Küpper et al. , 2015Mastrobuono-Battisti et al. 2013;Sanders & Binney 2013a,b;Bovy 2014;Fardal et al. 2015;Thomas et al. 2016;Ibata et al. 2020).
Also, we do not take into account the host system of the progenitor system that could affect the morphology of the GD-1 stream (Malhan et al. 2022;Qian et al. 2022). If we were to reproduce the density profile of the GD-1 stream as a function of φ 1 , the effects mentioned above are important and the only way to faithfully take these effects into account is to run N -body simulations or particle spray simulations (see Appendix D). However, because we are interested in the probability that the gaps were formed by the GCs, our approach is good enough for our purpose.
In this paper, we assumed that the MW potential is static and axisymmetric. These assumptions are simplistic, given that the MW is growing in time due to mass accretion (Buist & Helmi 2015), that the MW has a rotating bar (Hattori et al. 2016;Price-Whelan et al. 2016), and that the Large Magellanic Cloud has been perturbing the MW (Besla et al. 2010;Erkal et al. 2019;Garavito-Camargo et al. 2019;Koposov et al. 2019;Conroy et al. 2021;Petersen & Peñarrubia 2021;Shipp et al. 2021). These effects can alter the values of P gap-forming for each GC. However, since we did not tune the MW potential to maximize or minimize P gap-forming (instead, we just adopted one of the widely-used MW model potentials in our main analysis), our estimate of the probability that the GCs formed all three gaps in the GD-1 stream, 1.7 × 10 −5 (see equation (12)), is probably not too far from reality. Indeed, our additional analysis in Section 5.3, in which we varied the potential, supports this view. Thus, even though our simulation neglects some important physics, we believe our main conclusion is robust: the probability that GCs are responsible for all three gaps in the GD-1 stream is extremely low.
6. CONCLUSION In this paper, we estimated the probability that Galactic GCs can form a gap in the GD-1 stream by using test-particle simulations. The summary of this paper is as follows.
• There is a moderate trend that more massive GCs tend to have a larger P gap-forming (Fig. 3). However, the relative distance and velocity at their closest approach to the GD-1 stream are much more critical factors than the mass of the GCs.
• As shown in Figs. 6 and 7, our perturbed models can capture some of the morphological properties of the observed GD-1 stream, such as the length, widths, and strength of the gaps.
• The probability that all three gaps in the GD-1 stream were formed by the GCs is extremely low. In our fiducial model with T = 6 Gyr, this probability is P = 1.7 × 10 −5 (Section 5.1). This probability decreases if we adopt smaller T (see equation (14)). Assuming T ≤ 3 Gyr results in P = 0, which explains the result of Bonaca et al. (2019) who assumed T = 1 Gyr.
• The expected number of gaps in the GD-1 stream due to the flyby of GCs is N gap (GCs) = 0.057 in our fiducial model (equation (13)). This number is smaller than that due to the flyby of subhalos (N gap (subhalos) = 0.6) by a factor of 10 ).
• Given (i) that the probability that all three gaps in the GD-1 stream are formed by the GCs is extremely low, and (ii) that the retrograde orbit of the GD-1 stream makes other baryonic perturbers (e.g., spiral arms, the Galactic bar, or giant molecular clouds) even less likely to form the gaps, at least one of the gaps in the GD-1 stream is probably formed by the dark matter subhalos (Carlberg 2009(Carlberg , 2016Yoon et al. 2011;Erkal & Belokurov 2015a,b;Erkal et al. 2016;Bonaca et al. 2019;Banik et al. 2021).
• To sophisticate our analysis, we need to run Nbody simulations of the encounters of the GD-1 stream and GCs by also including the effect from the Large Magellanic Cloud (Erkal et al. 2019;Koposov et al. 2019;Shipp et al. 2021), or from the host halo of the GD-1's progenitor (Malhan et al. 2021;Qian et al. 2022).
The authors thank the referee for thorough reading and constructive comments that improved the original manuscript. DY and KH thank NAOJ for financial aid during the 2021 summer student program. We thank Junichi Baba for sharing his N -body code that gave insights into our work. KH thanks lecturers of N -body winter school 2021 held by NAOJ for stimulating lectures. KH is supported by JSPS KAKENHI Grant Numbers JP21K13965 and JP21H00053.
Numerical computations were in part carried out on GRAPE system at Center for Computational Astrophysics, National Astronomical Observatory of Japan. This research was supported in part through computational resources and services provided by Advanced Research Computing (ARC), a division of Information and Technology Services (ITS) at the University of Michigan, Ann Arbor. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/ consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration.
Facility: Gaia, LAMOST, SDSS/SEGUE Software: Agama(Vasiliev2019), corner.py(Foreman-APPENDIX A. COORDINATE SYSTEM We adopt a right-handed Galactocentric Cartesian coordinate system (x, y, z), which is the same as in Hattori et al. (2021). The position of the Sun is assumed to be x = (x , y , z ) = (−R 0 , 0, 0), with R 0 = 8.178 kpc (Gravity Collaboration et al. 2019). The velocity of the Sun with respect to the Galactic rest frame is assumed to be v = (v x, , v y, , v z, ) = (11.10, 247.30, 7.25) km s −1 (Reid & Brunthaler 2004;Schönrich et al. 2010). Following Koposov et al. (2010), we define the GD-1 coordinate (φ 1 , φ 2 ) such that The coordinate φ 1 is aligned with the track of the GD-1 stream, while the coordinate φ 2 is perpendicular to the stream. Due to the historical reason, φ 1 decreases along the direction of the GD-1's motion.
B. ASSESSMENT OF GD-1-LIKE GAPS As mentioned in Section 3.6, we judge whether a given perturbed stream model has a GD-1-like gap by comparing the perturbed model with the unperturbed model. Fig. 8 illustrates how we estimate the strength of the gap (under-dense region) in a model stream. In this example, we use the show-case model FSR1758 731. This perturbed model has a GD-1-like gap at φ 1,gap = −35 deg (which satisfies condition (B)). The density ratio at φ 1 = φ 1,gap is ρ perturbed /ρ unperturbed = 0.35, which is below our threshold of 0.8 (condition (A)). When T = 1 Gyr and 2 Gyr are assumed, no GCs form a gap in the GD-1 stream. When T = 3 Gyr and 4 Gyr are assumed, we found that two and five GCs form a gap, respectively, as summarized in Table 3.

D. COMMENTS ON THE PARTICLE SPRAY METHOD
In this paper, we release stream stars at once at t = −T , from the center of the progenitor. In reality, however, stream stars continuously escape from the progenitor from the neighborhood of the inner and outer Lagrange points. In order to check how our simple prescription affects the results, we additionally run eight simulations with the particle spray method, mostly following the prescription in Fardal et al. (2015) (and partially following Küpper et al. 2010) but with a constant stripping rate at −6 Gyr < t < t end with t end = (−5, −4, −3, −2) Gyr. In what follows, we focus on the show-case model NGC2808 698, in which the GD-1 stream and NGC2808 experience a close encounter at t = −4.154 Gyr. For each t end , we run two simulations where the mass of NGC2808 is set to be either 0 (unperturbed model) or 8.64 × 10 5 M (perturbed model). The result of these simulations is presented in Fig. 9.
The left-hand panels in Fig. 9 indicate that the unperturbed streams show a weak tendency such that the stream is more concentrated in φ 1 -direction if the disruption is completed more recently (smaller |t end |). This result can be intuitively understood: If a progenitor of the stream experience a more prolonged disruption history (smaller |t end |), we expect more stars near the progenitor at the current epoch (t = 0).
In the right-hand panels in Fig. 9, there is a visible gap at φ 1 −30 deg for simulations with t end = −5, −4, and −3 Gyr, while the gap is almost invisible for the simulation with t end = −2 Gyr. This result can be understood as follows. If t end = −2 Gyr, a large fraction of stream stars is still bound to the progenitor at t = −4.154 Gyr, when the GD-1 stream and NGC2808 experience a close encounter. Thus, the number of stars strongly affected by this perturbation is small. In addition, some fraction of stars that leave the progenitor after t = −4.154 Gyr can reach the gap region and 'fill' the gap as the stream evolves. These two effects make the contrast of the gap weaker for simulations with smaller |t end |.
The results of the particle spray method indicate that the gap is clearer for a larger value of |t end |, or more bursty stripping history. In this regard, our fiducial simulations -which corresponds to a bursty stripping or instantaneous disruption -tend to show clearer gaps. Therefore, our estimates of the gap-forming probability or the expected number of gaps in the main part of this paper may be an upper limit. Given that we find it extremely rare for GCs to form three gaps in the GD-1 stream by using a simplified simulations, our additional, more realistic simulations in this Appendix further support our conclusion.
Lastly, we comment on the difference between the particle spray simulations and our fiducial simulations. The streams generated from the particle spray method are more concentrated in φ 1 -direction than the streams gen-erated in our main analysis of this paper. This finding indicates that our simplified prescription in the fiducial analysis should not be used to interpret the global density profile ρ(φ 1 ) of the GD-1 stream. (However, we stress that our fiducial models are good enough to understand the gap-forming probability.) Figure 9. Additional simulations of the GD-1 stream model NGC2808 698 in which we used the particle spray method. In each simulation, we release the stars from the inner and outer Lagrange points with a constant stripping rate at −6 Gyr < t < t end . From top to bottom, we show the results with t end / Gyr = −5, −4, −3, and −2. The left column corresponds to the simulations with no perturbation. The right column corresponds to the simulations with perturbation from NGC2808. The gap becomes clearer if the disruption ends at an earlier epoch.