On N-body simulations of globular cluster stream
MMNRAS , 1–6 (0000) Preprint 29 January 2021 Compiled using MNRAS L A TEX style file v3.0 On 𝑁 -body simulations of globular cluster streams Nilanjan Banik ★ & Jo Bovy Mitchell Institute for Fundamental Physics and Astronomy, Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON, M5S 3H4, Canada
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Stellar tidal streams are sensitive tracers of the properties of the gravitational potential in which they orbit and detailedobservations of their density structure can be used to place stringent constraints on fluctuations in the potential caused by,e.g., the expected populations of dark matter subhalos in the standard cold dark matter paradigm (CDM). Simulations of theevolution of stellar streams in live 𝑁 -body halos without low-mass dark-matter subhalos, however, indicate that streams exhibitsignificant perturbations on small scales even in the absence of substructure. Here we demonstrate, using high-resolution 𝑁 -bodysimulations combined with sophisticated semi-analytic and simple analytic models, that the mass resolutions of 10 –10 M (cid:12) commonly used to perform such simulations cause spurious stream density variations with a similar magnitude on large scalesas those expected from a CDM-like subhalo population and an order of magnitude larger on small, yet observable, scales. Weestimate that mass resolutions of ≈
100 M (cid:12) ( ≈ (cid:12) ) are necessary for spurious, numerical density variations to be well belowthe CDM subhalo expectation on large (small) scales. That streams are sensitive to a simulation’s particle mass down to suchsmall masses indicates that streams are sensitive to dark matter clustering down to these low masses if a significant fraction ofthe dark matter is clustered or concentrated in this way, for example, in MACHO models with masses of 10–100 M (cid:12) . Key words:
Cosmology: dark matter — Galaxy: evolution — Galaxy: halo — Galaxy: kinematics and dynamics — Galaxy:structure
The standard cold dark matter (CDM) picture of structure formationpredicts that a Milky-Way-size galaxy halo should host numerousdark substructures (“subhalos”) with a range of masses down tomany of orders of magnitudes below the typical dwarf galaxy mass(Diemand et al. 2008; Springel et al. 2008). Detecting these substruc-tures can not only provide tell-tale evidence for the existence of darkmatter but will also enable us to put stringent constraints on its parti-cle nature. Unfortunately, being low in mass, these substructures arenot able to initiate star formation and hence are undetectable directlyby telescopes observing in the electromagnetic spectrum.Flux perturbations in strong gravitational lensing systems due tothese substructures provides an indirect way of detecting and analyz-ing these low mass substructures (Dalal & Kochanek 2002; Vegetti& Koopmans 2009; Despali & Vegetti 2017; Gilman et al. 2020).A powerful complementary approach is to study stellar density per-turbations along stellar streams left as a result of gravitational en-counters with these dark subhalos. Stellar streams emerging fromtidally disrupting globular clusters in our Galaxy have a largely one-dimensional structure and in the absence of external perturbationshave a fairly uniform stellar density along its length (Johnston et al.1999; Sanders & Binney 2013; Bovy 2014). Owing to the very lowvelocity dispersion among its member stars, these streams are dy-namically cold and hence are extremely sensitive to perturbations ★ E-mail: [email protected] in the underlying gravitational potential. Close flybys of subhaloscan therefore leave visible imprints in the stellar density along thestreams in the form of gaps (Ibata et al. 2002; Johnston et al. 2002;Siegal-Gaskins & Valluri 2008; Carlberg 2009). While these gaps cancan be individually analyzed to infer the properties of the perturbingsubhalo (Yoon et al. 2011; Carlberg 2012, 2013; Erkal & Belokurov2015a,b; Sanders et al. 2016), powerful statistical techniques basedon the power spectrum of the full stream density (Bovy et al. 2017)have been applied to observed GD-1 and Pal 5 stream data to con-strain the abundance of subhalos in the mass range 10 − M (cid:12) and constrain the particle mass of thermal dark matter (Banik et al.2019, 2020).A key step in effectively using streams as probes for dark mattersubstructure is to identify sources of stream density perturbationsother than subhalo encounters and modeling their effects in the anal-ysis. Depending on the orbit, stellar streams can be significantlyperturbed by the baryonic structures in our Galaxy such as the bar(Erkal et al. 2017; Pearson et al. 2017; Banik & Bovy 2019), the spiralarms (Banik & Bovy 2019) and the Giant Molecular Clouds (GMCs)(Amorisco et al. 2016; Banik & Bovy 2019). Another source ofstream density variations are the epicyclic overdensities which arisesince stars are tidally ejected in bursts near the pericentric passageof the progenitor globular cluster (Kuepper et al. 2010, 2012). Suchoverdensities along the stream due to episodic tidal stripping werestudied in (Sanders et al. 2016; Bovy et al. 2017) and were found toquickly disperse out as the ejected stars mixed due to their velocitydispersion. Stream regions near the progenitor that are dynamically © a r X i v : . [ a s t r o - ph . GA ] J a n Banik & Bovy young may still have such overdensities since the stars did not havesufficient time to mix, therefore such regions need to be judiciouslyexcluded in inferring dark matter properties.Recently, Carlberg (2018) and Carlberg (2020) investigated thedynamical evolution of globular clusters and their tidally disruptedstreams in cosmological 𝑁 -body simulations by placing the progeni-tor clusters on near circular orbits inside reconstituted subhalos fromthe Via-Lactea II simulation (Madau et al. 2008) and evolving them.These subhalos along with their globular clusters and their tidallydisrupted streams were eventually accreted onto the host halo andthe final stream structures had a wide range of density variations. Inorder to test whether these density variations were caused by grav-itational encounters with lower mass subhalos ( (cid:46) × M (cid:12) ), asimilar simulation without these lower mass subhalos was run andthe resulting stream densities still had similar variations. Based onthis finding, it was concluded that the stream density variations werenot caused by the low mass subhalo impacts.However, while Carlberg (2018) was careful to use a mass reso-lution of 2 × M (cid:12) such that heating from the massive 𝑁 -bodyparticles is insignificant, it is unclear whether coherent perturbationsto the kinematics of simulated tidal streams by the 𝑁 -body particlesare significant. Coherent perturbations are what is relevant for deter-mining whether stream density variations from numerical effects canbe mistaken for true subhalo-induced variations, because the latterproduce only coherent fluctuations without significant heating. Be-cause 𝑁 -body simulations are useful for understanding the evolutionof tidal streams in the full cosmological context, in this paper, weinvestigate whether significant stream density variations are createdin simulated globular-cluster streams evolved within a Milky Waylike live halo of different mass resolutions (10 M (cid:12) , similar to Carl-berg 2020, and 10 M (cid:12) ). We show that in the absence of subhalos,numerically-induced stream density fluctuations at these mass reso-lutions exceed the expected subhalo signal by an order of magnitudeand we derive the necessary mass resolution for numerical effects tobe well below the expected CDM subhalo signal. We run three different 𝑁 -body simulations of the dynamical evolutionof a globular cluster with an orbit that is similar to that of the GD-1 stream (Grillmair & Dionatos 2006) in the Milky Way, becauseGD-1 is currently the best stream to study subhalo-induced densityperturbations (Banik et al. 2019). We let the cluster and stream evolvefor 4 Gyr, the approximate age of the GD-1 stream. The progenitorcluster is modeled by an isothermal King profile (King 1966) withtotal mass 7360 M (cid:12) , a half-mass radius of 20 pc, and a dimensionlesscentral potential depth 𝑊 is set to 7. We use LIMEPY (Gieles &Zocchi 2018) to generate the initial condition of the cluster with100,000 equal-mass particles; because the GD-1 progenitor and itsstars are likely ≈
12 Gyr old, the effects of stellar evolution andrelaxation are small, allowing us to model the system as a collisionlesssystem (Webb & Bovy 2019). The velocity dispersion of the initialcluster is 0 .
46 km s − . The host potential consists of a Milky-Way-like halo and disk. The halo is modeled as an NFW profile (Navarroet al. 1997) using the best fit parameters from (McMillan 2017)(Table 3), i.e. scale length 𝑟 ℎ = . 𝜌 = . × M (cid:12) kpc − . We generate two live Milky Way likehalo initial conditions using GalactICS (Kuijken & Dubinski 1995;Widrow & Dubinski 2005; Widrow et al. 2008; Deg et al. 2019) withparticle mass resolution of 10 M (cid:12) and 10 M (cid:12) . The left panel ofFigure 1 shows that the rotation curve of the live halo is consistent R (kpc) v c ( k m / s ) Live haloAnalytic halo R (kpc) Interpolated haloLive halo Figure 1.
Rotation curves of our live, analytic, and smooth halos.
Left panel :the solid red curve shows the rotation curve of the live 𝑁 -body halo, whichis initialized as an NFW profile and has a mass resolution of 10 M (cid:12) , andthe dashed black curve shows the analytic rotation curve with the same NFWprofile parameters. Right panel : the rotation curve of the same live halois again shown as a solid red curve and the dashed black curve displays therotation curve of the static, smooth interpolated representation of the same livehalo as described in the text. All of the rotation curves agree well, indicatingthat the overall mass distribution of all of the models that we simulate issimilar. with that of the analytic NFW halo with the same parameters. Thedisk is modeled as a Miyamoto-Nagai profile with a total mass of6 . × M (cid:12) , a radial scale length of 3 kpc, and a scale height of 0.28kpc (Bovy 2015). Because we are mainly interested in the studying theeffects of mass resolution of the halo on the globular cluster stream,we evaluate the disk potential analytically as an external force in the 𝑁 -body simulations. To study the evolution of the globular clusterstream in a smooth halo similar to the 𝑁 -body one, we determinethe radial acceleration 𝑎 𝑟 ( 𝑟 ) by computing − 𝐺 𝑀 ( < 𝑟 )/ 𝑟 on auniform radial grid using the positions of the live halo particles in the10 M (cid:12) resolution case and interpolating them to create a spherical,static, and smooth representation of the live halo. The right panelof Figure 1 compares the rotation curve of the live halo with thatfrom the interpolated smooth halo within 50 kpc of the host centerand shows good agreement. We limit our interpolation to 50 kpc,because the stream orbit is confined well within this radial range. Tosummarize: we run three main simulations: a live halo with massresolution 10 M (cid:12) , a live halo with 10 M (cid:12) resolution, and a statichalo created from smoothing out the 10 M (cid:12) resolution live halothrough interpolation.We obtain a GD-1 like orbit for the cluster for each of the threesimulations that we run by placing our cluster at the assumed cur-rent phase space coordinates of the GD-1 progenitor following Webb& Bovy (2019), flipping the sign of the velocity, and backwards-evolving for 4 Gyr using the 𝑁 -body code GIZMO (Hopkins 2015)which is based on the GADGET code (Springel 2005). In all simu-lations, we use a force softening of 2 pc for the cluster particles andin the live-halo simulations we use a softening of 200 pc and 20 pcfor the halo particles in the 10 M (cid:12) and 10 M (cid:12) resolution cases,respectively. At the end of the backwards-evolution, we determinethe phase space coordinate of the cluster’s center of density using clustertools and we then place the same initial cluster at thispoint, flip the sign of the velocity, and forward-evolve it for 4 Gyr.At the end of the backward-forward simulation, the resulting clustercenter was found to match the present day GD-1 progenitor’s phasespace location that we started with very well and a realistic tidalstream of escaped cluster particles forms in each simulation. The Available at https://github.com/webbjj/clustertools .MNRAS000
Left panel :the solid red curve shows the rotation curve of the live 𝑁 -body halo, whichis initialized as an NFW profile and has a mass resolution of 10 M (cid:12) , andthe dashed black curve shows the analytic rotation curve with the same NFWprofile parameters. Right panel : the rotation curve of the same live halois again shown as a solid red curve and the dashed black curve displays therotation curve of the static, smooth interpolated representation of the same livehalo as described in the text. All of the rotation curves agree well, indicatingthat the overall mass distribution of all of the models that we simulate issimilar. with that of the analytic NFW halo with the same parameters. Thedisk is modeled as a Miyamoto-Nagai profile with a total mass of6 . × M (cid:12) , a radial scale length of 3 kpc, and a scale height of 0.28kpc (Bovy 2015). Because we are mainly interested in the studying theeffects of mass resolution of the halo on the globular cluster stream,we evaluate the disk potential analytically as an external force in the 𝑁 -body simulations. To study the evolution of the globular clusterstream in a smooth halo similar to the 𝑁 -body one, we determinethe radial acceleration 𝑎 𝑟 ( 𝑟 ) by computing − 𝐺 𝑀 ( < 𝑟 )/ 𝑟 on auniform radial grid using the positions of the live halo particles in the10 M (cid:12) resolution case and interpolating them to create a spherical,static, and smooth representation of the live halo. The right panelof Figure 1 compares the rotation curve of the live halo with thatfrom the interpolated smooth halo within 50 kpc of the host centerand shows good agreement. We limit our interpolation to 50 kpc,because the stream orbit is confined well within this radial range. Tosummarize: we run three main simulations: a live halo with massresolution 10 M (cid:12) , a live halo with 10 M (cid:12) resolution, and a statichalo created from smoothing out the 10 M (cid:12) resolution live halothrough interpolation.We obtain a GD-1 like orbit for the cluster for each of the threesimulations that we run by placing our cluster at the assumed cur-rent phase space coordinates of the GD-1 progenitor following Webb& Bovy (2019), flipping the sign of the velocity, and backwards-evolving for 4 Gyr using the 𝑁 -body code GIZMO (Hopkins 2015)which is based on the GADGET code (Springel 2005). In all simu-lations, we use a force softening of 2 pc for the cluster particles andin the live-halo simulations we use a softening of 200 pc and 20 pcfor the halo particles in the 10 M (cid:12) and 10 M (cid:12) resolution cases,respectively. At the end of the backwards-evolution, we determinethe phase space coordinate of the cluster’s center of density using clustertools and we then place the same initial cluster at thispoint, flip the sign of the velocity, and forward-evolve it for 4 Gyr.At the end of the backward-forward simulation, the resulting clustercenter was found to match the present day GD-1 progenitor’s phasespace location that we started with very well and a realistic tidalstream of escaped cluster particles forms in each simulation. The Available at https://github.com/webbjj/clustertools .MNRAS000 , 1–6 (0000) n 𝑁 -body simulations of globular cluster streams [] (a) No subhalos, M (cid:12) resolution N -body[] (b) No subhalos, M (cid:12) resolution N -body φ ( d e g ) [] (c) No subhalos, smoothed N -body − − − −
25 0 25 50 75 100 φ (deg) (d) CDM subhalos Figure 2.
The simulated streams after 4 Gyr of evolution represented in the ( 𝜙 , 𝜙 ) angular sky coordinate frame centered at the progenitor location. Panel (a): 𝑁 -body run with live halo of mass resolution of 10 M (cid:12) , panel (b): 𝑁 -body run with live halo of mass resolution of 10 M (cid:12) , panel (c): static,smoothed 𝑁 -body halo, panel (d): simulated stream in the streampepperdf framework that was impacted by the CDM abundance of subhalos in the massrange 10 − − M (cid:12) . For the 𝑁 -body runs, the progenitor can still be seenat 𝜙 ∼
0. In both panels (a) and (b), significant density variations can beseen along the stream. In panel (c), epicylic overdensities can be clearly seenwithin 20 ◦ around the progenitor. Overall, the streams evolved in the live 𝑁 -body realization of a smooth halo show significantly more density variationsthan in the smooth or CDM-like halos. mean Galactocentric radius of the orbit is ∼
17 kpc, with a meanorbital eccentricity of 0.2, and peri- and apogalacticon distances of ∼
13 kpc and ∼
20 kpc, respectively, similar to the observed GD-1stream (Webb & Bovy 2019). Because the initial condition of thehalo was generated in the absence of a disk, placing the analytic diskat the center of the halo breaks its state of equilibrium. Therefore,we implement conditions on the analytic disk such that the halo par-ticles (particle type 1) do not feel the gravitational force from thedisk, while the cluster particles (particle type 4) feel the combinedforce from the halo and the disk. We have also ignored adiabaticcontraction of the halo which should not have any effect at the smallscales of stream density variations. 𝑁 -BODY, SMOOTH,AND CDM HALOS The resulting streams in all of our simulations are shown in panels(a),(b) and (c) in Figure 2 where ( 𝜙 , 𝜙 ) are coordinates along andtransverse to the stream centered at the progenitor. The ( 𝜙 , 𝜙 ) coordinate system is a rotation of the sky coordinate system that wasdetermined by hand to align the stream to have constant 𝜙 suchthat 𝜙 can be used as the along-stream angular coordinate. In bothlive halo cases the resulting streams can be seen to have acquiredsubstantial density variations that are clearly absent in the smoothedhalo case. For reference, we show in panel (d) one realization of asimilar stream that was impacted by a CDM abundance of subhalosin the mass range 10 − M (cid:12) . This last case is generated using the [] (a) No subhalos, M (cid:12) resolution N -bodyTrailing Leading (b) No subhalos, M (cid:12) resolution N -body (c) No subhalos, smoothed N -body N o r m a li z e d D e n s i t y − − − − − − φ (deg) (d) CDM subhalos
20 25 30 35 40 45 φ (deg) Figure 3.
Normalized stream density of the four cases shown in Figure 2.This is computed by first binning the particles along the stream in 0 . ◦ binsin 𝜙 and then dividing by a 3 rd order smoothing polynomial fit. The rightcolumns show the leading arm and the left columns show the trailing armafter excluding 20 ◦ around the progenitor and considering 30 ◦ and 25 ◦ alongthe leading and trailing arm respectively. The red error bars are the shot noisein each bin. − − − − − − φ (deg) . . . . . . . σ v ( k m / s ) Trailing
20 25 30 35 40 45 φ (deg)LeadingCDM subhalosNo subhalos, smoothed N -bodyNo subhalos, M (cid:12) resolution N -bodyNo subhalos, M (cid:12) resolution N -body Figure 4.
Velocity dispersion along the stream in the four different casesshown in Figure 2. There is no substantial overall heating along the stream inany of the simulated cases. stream-subhalo interaction framework streampepperdf , which isbased on the galpy code (Bovy 2015). The CDM subhalo abundanceand their corresponding sizes were obtained following the fittingfunctions from Erkal & Belokurov (2015b) that were based on theAquarius simulations (Springel et al. 2008). The subhalos are set toimpact the stream following the same steps as in Bovy et al. (2017)and we refer the reader to that paper for full details.To quantify the amount of density variations in each case, wewill compute the power spectrum of the normalized stream densityfollowing Bovy et al. (2017). Variations in the stream density dueto epicyclic pile ups can be most clearly seen within ∼ ◦ around Available at https://github.com/jobovy/streamgap-pepper .MNRAS , 1–6 (0000)
Banik & Bovy /k φ (deg) . . p P δδ ( k φ ) Trailing M res = 10 M (cid:12) N -body M res = 10 M (cid:12) N -bodySmoothed N -bodyCDM /k φ (deg) Leading
Figure 5.
Density power spectrum of the leading and trailing arm of streams simulated in live 𝑁 -body halos with 10 M (cid:12) (red curve) and 10 M (cid:12) resolution.The green curve shows the power spectrum of the stream evolved in the static, smoothed 𝑁 -body run. The black curve and the gray shaded region shows themedian and 2 𝜎 dispersion of density power due to CDM subhalo impacts. At the largest scales the density power in the live halo cases are comparable to thatdue to impacts by CDM subhalos, while for the smoothed halo case the density power is substantially lower and consistent with noise. On small scales, thedensity power in the live halo cases exceed that expected from a CDM-like population of subhalos by an order of magnitude or more. the progenitor in the smoothed halo case, therefore we cut out thatregion from each case to exclude their contribution to the densityvariations. Furthermore, we restrict our analysis to 27 ◦ and 25 ◦ alongthe leading and trailing arms respectively for each case. This way weare comparing the density variations over the same angular rangealong the stream for the different cases. The densities are normalizedby dividing out by a 3 rd order smoothing polynomial fit following(Bovy et al. 2017). The resulting normalized stream densities areshown in Figure 3. As expected from Figure 2, the stream density inFigure 3 in the smoothed halo case is mostly flat, while the live halocases display wide-scale density variations. The velocity dispersionalong the stream length for the different cases are shown in Figure4, which demonstrates that heating of the less massive star particlesof the cluster by the more massive dark matter particles of the halois not significant, consistent with the analysis presented in (Carlberg2018).The power spectrum of the stream density, computed using thesame technique from (Bovy et al. 2017), from the live and smoothedhalo 𝑁 -body runs are shown in Figure 5 by the red solid curve(10 M (cid:12) resolution), the blue solid curve (10 M (cid:12) resolution), andthe green solid curve (smoothed static halo). Overplotted as the blacksolid line and gray shaded region are the median and 2 𝜎 dispersionin the stream density power respectively, of 1000 realizations of thestream impacted by a CDM abundance of subhalos. From this figureit is clear that at the largest angular scales, the density power accruedby the stream through the interactions with the halo particles in the 𝑁 -body runs is comparable to that due to the CDM subhalo impacts.The power in the smoothed 𝑁 -body halo case has comparatively verylow power that is consistent with Poisson noise, which is ≈ .
04 forboth arms. The power in both the live halo cases is at all scales muchhigher than the noise. This shows that even in the absence of subhaloimpacts, interactions with the dark matter halo particles of resolution10 M (cid:12) or 10 M (cid:12) can substantially perturb the stream density. Next, to conclusively demonstrate that these stream density vari-ations are really due to interactions with the 𝑁 -body particles, weevolve the same stream in the smoothed halo potential while impact-ing it with a population of subhalos similar in abundance and struc-ture as the 𝑁 -body particles. We want to check whether the resultingstream density power spectrum is comparable to those obtained inthe live halo cases. To do this, we again use the fast stream-subhalointeraction simulator streampepperdf and model the subhalos asPlummer spheres with mass equal to the 𝑁 -body particle mass andscale length equal to the Plummer equivalent softening of the 𝑁 -bodyparticles. The latter in GIZMO is 1 / . ∼
17 kpc in our case. The numberdensity of 𝑁 -body particles in the 10 M (cid:12) resolution case is around2513 times the CDM predicted subhalo number density in the massrange 10 . − . M (cid:12) and for the 10 M (cid:12) resolution case this factor 𝑛 𝑁 − body / 𝑛 subhalo is around 3168 times the CDM predicted subhalonumber density in the mass range 10 . − . M (cid:12) . Simulatingall encounters from such a huge number of subhalos is computa-tionally very challenging. Therefore we use an approximate methodin which we evolve the stream with 10 × the CDM predicted sub-halo number density and scale the resulting density power spectrumby 𝑛 𝑁 − body / 𝑛 subhalo /
10, because uncorrelated subhalo encountersshould contribute additively to the power spectrum (note that whatwe plot in Figure 5 is the square root of the power spectrum). Weran 1000 simulations each for the 10 M (cid:12) and 10 M (cid:12) resolutioncases and the results are shown in Figure 6. The red and blue dashedcurves show the median power in the 10 M (cid:12) and 10 M (cid:12) resolutioncases, respectively. The 2 𝜎 dispersion of power in the correspondingcases is shown by the shaded regions. The density power from thelive halo cases (solid curves) are consistent with the correspondingpredicted density powers within their dispersion confirming that the MNRAS000
10, because uncorrelated subhalo encountersshould contribute additively to the power spectrum (note that whatwe plot in Figure 5 is the square root of the power spectrum). Weran 1000 simulations each for the 10 M (cid:12) and 10 M (cid:12) resolutioncases and the results are shown in Figure 6. The red and blue dashedcurves show the median power in the 10 M (cid:12) and 10 M (cid:12) resolutioncases, respectively. The 2 𝜎 dispersion of power in the correspondingcases is shown by the shaded regions. The density power from thelive halo cases (solid curves) are consistent with the correspondingpredicted density powers within their dispersion confirming that the MNRAS000 , 1–6 (0000) n 𝑁 -body simulations of globular cluster streams /k φ (deg) . . p P δδ ( k φ ) Trailing /k φ (deg) Leading M res = 10 M (cid:12) N -body M res = 10 M (cid:12) N -body M res = 10 M (cid:12) N -body equiv. M res = 10 M (cid:12) N -body equiv. Figure 6.
Power spectrum of the leading and trailing arm in the live 𝑁 -body runs (solid curves) compared with the case where the stream was impacted by aperturber population that has a similar abundance as the 𝑁 -body particles and for which the perturbers have the same mass and similar structure as the softened 𝑁 -body particles. For the latter case, the median (dashed curve) and 2 𝜎 dispersion (shaded region) are shown. The red color is used for the 10 M (cid:12) casewhile the blue color for the 10 M (cid:12) case. The good agreement between the 𝑁 -body results and the perturber-population prediction demonstrates that the massresolution is the dominant reason for the large density variations in the 𝑁 -body streams in Figure 2. density variations in the live 𝑁 -body simulations are artefacts of thedark matter particle resolution. It can also be seen that for both thesecases as well as for the live halo cases, there is much more powerat small scales compared to the CDM case, which is expected sinceas shown in Bovy et al. (2017), low mass (such as ∼ M (cid:12) ) per-turbers result in small scale density fluctuations. That said, we alsofind enormous power at large scales along the stream in these caseswhich seems counter-intuitive. This is because density fluctuationsinflicted on streams grow linearly with time (Erkal & Belokurov2015a) and so over the course of its evolution the stream densitystructure is a combination of older large scale fluctuations and recentsmall scale fluctuations. We have investigated the effects of dark-matter particle mass res-olution in 𝑁 -body simulations on the evolution of globular clusterstreams. Comparing 𝑁 -body simulations of the evolution of a glob-ular cluster and its tidal stream in live Milky Way like halos withresolutions of 10 M (cid:12) and 10 M (cid:12) with the same stream evolved ina static, smoothed dark matter halo demonstrates that, despite the ab-sence of significant overall heating along the stream, in the live halosimulations the final stream density displays significant variations onall angular scales along the stream. In particular, the density powerat the largest scales is comparable to that expected due to impacts bya CDM-like population of subhalos in the mass range 10 –10 M (cid:12) .We further verify that the stream density variations are due to inter-actions with the massive 𝑁 -body particles by evolving the stream inthe smoothed halo while impacting it with a population of subha-los similar in abundance and structure as the 𝑁 -body particles. Wefound that the resulting density power is in agreement with that ofthe live halo cases. This conclusively shows that 𝑁 -body resolutions of 10 M (cid:12) and 10 M (cid:12) can severely perturb globular cluster streamdensities and that these mass resolutions—while high by present-daystandards—are inadequate for studying the expected density varia-tions along tidal streams in 𝑁 -body simulations. In particular, theresults presented here demonstrate that the density variations foundin the cosmological simulations of Carlberg (2020) in the absence ofsubhalo impacts are likely a resolution artifact owing to using a darkmatter resolution of 2 × M (cid:12) .While we have used 𝑁 -body simulations to determine the expectedpower spectrum of stream density variations in the presence of apopulation of massive 𝑁 -body particles, when the entire halo ismade up off 𝑁 = 𝑀 halo / 𝑀 res particles with mass 𝑀 res , the number ofimpacts on a stream is high enough that we can estimate the inducedpower using simple statistical calculations. Following Appendix B ofDalal et al. (2020), who work out the expected stream density powerspectrum for a population of CDM subhalos, we can determine theapproximate scaling of the power spectrum with mass. As shown byDalal et al. (2020), for perturbers of a single mass 𝑀 res with a spatialdensity ¯ 𝑛 , the expected power spectrum is 𝑃 𝛿 𝛿 ∝ 𝑀 ¯ 𝑛 . When theentire halo consists of substructure of mass 𝑚 , then ¯ 𝑛 = 𝜌 / 𝑀 res ,where 𝜌 is the dark matter density and 𝑃 𝛿 𝛿 ∝ 𝑀 res (where we dropthe 𝜌 dependence because it does not change with resolution). Thisscaling holds down to angular scales 𝜃 ≈ 𝑅 / 𝑟 , where 𝑅 is the size ofthe perturber and 𝑟 is the Galactocentric radius (that is, 𝜃 is the angularsize of perturbers as seen from the Galactic center); below this scalethe power spectrum drops to zero quickly. Figures 5 and 6 show thaton large scales indeed approximately 𝑃 𝛿 𝛿 ∝ 𝑀 res . To be able tosimulate the expected large-scale power without numerical artefactsrequires at least an order of magnitude improvement in the 𝑁 -body √ 𝑃 𝛿 𝛿 or a decrease in 𝑁 -body particle mass of 100 (to 𝑀 res ≈
100 M (cid:12) ) and ideally smaller to avoid statistical upwards fluctuationsin the numerical power. Because the power spectrum of artificialfluctuations declines much less steeply than that predicted by CDM
MNRAS , 1–6 (0000)
Banik & Bovy subhalos, the resolution requirement becomes even more stringent atsmall scales; on the smallest scales at which the power is observablewith future data (few degree; Bovy et al. 2017) the numerical power √ 𝑃 𝛿 𝛿 is an order of magnitude larger than the predicted CDM powerand a mass resolution of 1 M (cid:12) would be necessary for artefacts tobe negligible. Such mass resolutions are higher than typically used,but are achievable with some current algorithms making efficient useof modern high-performance computing such as GPUs (e.g., Bédorfet al. 2014; Asano et al. 2020).Our results confirm that stellar streams in the Milky Way are ex-tremely sensitive to fluctuations in the gravitational potential. While,as we have shown, this puts stringent constraints on 𝑁 -body simula-tions of the evolution of such streams, it also presents an opportunityin that it further confirms that streams present one of the most sen-sitive ways to learn about non-standard dark matter models. In par-ticular, our results indicate that streams are sensitive to dark matterclustering on mass scales of ≈ (cid:12) if a significant fraction ofthe dark matter participates in this clustering. One such dark-mattercandidate is MACHOs (e.g., primordial black holes), where next-generation stream-density measurements on small scales by LSST atthe Vera Rubin Observatory (Ivezić et al. 2019) or the Roman SpaceTelescope (Spergel et al. 2015) may be able to constrain MACHOdark matter in the interesting ≈ (cid:12) range (Bird et al. 2016;Brandt 2016; Carr & Kühnel 2020). ACKNOWLEDGEMENTS
We thank Jeremy Webb and Nathan Deg for useful discussions. NBthanks Nathan Deg for help in setting up and running GalactICS. JBreceived financial support from NSERC (funding reference numberRGPIN-2020-04712), an Ontario Early Researcher Award (ER16-12-061), and from the Canada Research Chair program. Portionsof this research were conducted with high performance researchcomputing resources provided by Texas A&M University ( https://hprc.tamu.edu ). DATA AVAILABILITY
No new data were generated in support of this research.
REFERENCES
Amorisco N. C., Gòmez F. A., Vegetti S., White S. D. M., 2016, MNRAS,463, L17Asano T., Fujii M. S., Baba J., Bédorf J., Sellentin E., Portegies Zwart S.,2020, MNRAS, 499, 2416Banik N., Bovy J., 2019, MNRAS, 484, 2009Banik N., Bovy J., Bertone G., Erkal D., de Boer T. J. L., 2019, arXiv e-prints,p. arXiv:1911.02662Banik N., Bovy J., Bertone G., Erkal D., de Boer T. J. L., 2020, arXiv e-prints,p. arXiv:1911.02663Bédorf J., Gaburov E., Fujii M. S., Nitadori K., Ishiyama T., Porte-gies Zwart S., 2014, in Proceedings of the International Conferencefor High Performance Computing. pp 54–65 ( arXiv:1412.0659 ),doi:10.1109/SC.2014.10Bird S., Cholis I., Muñoz J. B., Ali-Haïmoud Y., Kamionkowski M., KovetzE. D., Raccanelli A., Riess A. G., 2016, Phys. Rev. Lett., 116, 201301Bovy J., 2014, ApJ, 795, 95Bovy J., 2015, ApJS, 216, 29Bovy J., Erkal D., Sanders J. L., 2017, MNRAS, 466, 628Brandt T. D., 2016, ApJ, 824, L31 Carlberg R. G., 2009, ApJ, 705, L223Carlberg R. G., 2012, ApJ, 748, 20Carlberg R., 2013, ApJ, 775, 90Carlberg R. G., 2018, ApJ, 861, 69Carlberg R. G., 2020, ApJ, 889, 107Carr B., Kühnel F., 2020, Annual Review of Nuclear and Particle Science,70, annurevDalal N., Kochanek C. S., 2002, ApJ, 572, 25Dalal N., Bovy J., Hui L., Li X., 2020, arXiv e-prints, p. arXiv:2011.13141Deg N., Widrow L. M., Randriamampandry T., Carignan C., 2019, MNRAS,486, 5391Despali G., Vegetti S., 2017, MNRAS, 469, 1997Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J.,2008, Nature, 454, 735Erkal D., Belokurov V., 2015a, MNRAS, 450, 1136Erkal D., Belokurov V., 2015b, MNRAS, 454, 3542Erkal D., Koposov S. E., Belokurov V., 2017, MNRAS, 470, 60Gieles M., Zocchi A., 2018, MNRAS, 474, 3997Gilman D., Birrer S., Nierenberg A., Treu T., Du X., Benson A., 2020,MNRAS, 491, 6077Grillmair C. J., Dionatos O., 2006, ApJ, 643, L17Hopkins P. F., 2015, MNRAS, 450, 53Ibata R. A., Lewis G. F., Irwin M. J., 2002, MNRAS, 332, 915Ivezić Ž., Kahn S. M., Tyson J. A., et al., 2019, ApJ, 873, 111Johnston K. V., Zhao H., Spergel D. N., Hernquist L., 1999, ApJ, 512, L109Johnston K. V., Spergel D. N., Haydn C., 2002, ApJ, 570, 656King I. R., 1966, AJ, 71, 64Kuepper A. H. W., Kroupa P., Baumgardt H., Heggie D. C., 2010, MNRAS,401, 105Kuepper A. H. W., Lane R. R., Heggie D. C., 2012, MNRAS, 420, 2700Kuijken K., Dubinski J., 1995, MNRAS, 277, 1341Madau P., Diemand J., Kuhlen M., 2008, ApJ, 679, 1260McMillan P. J., 2017, MNRAS, 465, 76Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Pearson S., Price-Whelan A. M., Johnston K. V., 2017, Nature Astronomy, 1,633Sanders J. L., Binney J., 2013, MNRAS, 433, 1826Sanders J. L., Bovy J., Erkal D., 2016, MNRAS, 457, 3817Siegal-Gaskins J. M., Valluri M., 2008, ApJ, 681, 40Spergel D., Gehrels N., Baltay C., et al., 2015, arXiv e-prints, p.arXiv:1503.03757Springel V., 2005, MNRAS, 364, 1105Springel V., et al., 2008, MNRAS, 391, 1685Vegetti S., Koopmans L. V. E., 2009, MNRAS, 392, 945Webb J. J., Bovy J., 2019, MNRAS, 485, 5929Widrow L. M., Dubinski J., 2005, ApJ, 631, 838Widrow L. M., Pym B., Dubinski J., 2008, ApJ, 679, 1239Yoon J. H., Johnston K. V., Hogg D. W., 2011, ApJ, 731, 58This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000