Generation of mock tidal streams
MMon. Not. R. Astron. Soc. , 1–20 (2014) Printed 2 July 2018 (MN L A TEX style file v2.2)
Generation of mock tidal streams
Mark A. Fardal ∗ , Martin D. Weinberg , Shuiyao Huang Dept. of Astronomy, University of Massachusetts, Amherst, MA, 01003, USA
Submitted to MNRAS 2 July 2018
ABSTRACT
In this paper we discuss a method for the generation of mock tidal streams. Usingan ensemble of simulations in an isochrone potential where the actions and frequen-cies are known, we derive an empirical recipe for the evolving satellite mass and thecorresponding mass loss rate, and the ejection conditions of the stream material. Theresulting stream can then be quickly generated either with direct orbital integration,or by using the action-angle formalism. The model naturally produces streaky featureswithin the stream. These are formed due to the radial oscillation of the progenitor andthe bursts of stars emitted near pericenter, rather than clumping at particular oscil-lation phases as sometimes suggested. When detectable, these streaky features are areliable diagnostic for the stream’s direction of motion and encode other informationon the progenitor and its orbit. We show several tests of the recipe in alternate poten-tials, including a case with a chaotic progenitor orbit which displays a marked effecton the width of the stream. Although the specific ejection recipe may need adjustingwhen elements such as the orbit or satellite density profile are changed significantly,our examples suggest that model tidal streams can be quickly and accurately generatedby models of this general type for use in Bayesian sampling.
Key words: galaxies: kinematics and dynamics – galaxies: interactions – galaxies:haloes – galaxies: star clusters
Tidal streams are interesting as potentially very high-precision probes of the mass and shapes of galactic poten-tials. They also mark the deaths of galactic satellites (whichextends the sample of dwarf galaxies we can study), serve asa link between intact galaxies and the smoother halo compo-nent, test the predicted population of dark subhalos withingalaxy halos, and provide insight into the formation of glob-ular cluster systems.The methods used for modeling tidal streams range indifficulty from simple orbit fits to use of N -body models. Theaccuracy of the analysis method should bear some relationto the quality of the observational data and the visible com-plexity of the stream under discussion. Some tidal streamsappear simply as a slightly broadened track, while othersappear to have multiple components (as in the case of theSagittarius stream, Belokurov et al. 2006 and Sohn et al.2014) or density variations in excess of random noise (as inthe case of Pal 5, Odenkirchen et al. (2003). Obtaining a fitto a given stream can be useful, but it is more informative toobtain confidence intervals for parameters of interest, whichusually involves Bayesian sampling techniques that requiremany stream-model comparisons. ∗ E-mail: [email protected]
Some recent papers have proposed stream analysismethods that avoid specific modeling of the tidal disruptionprocess (Price-Whelan et al. 2014; Sanderson et al. 2014).These methods work well when given high precision informa-tion on all six dimensions of phase space. But in most cases,and especially when some dimensions are unconstrained byobservations, it seems useful to have the analysis methodincorporate more information on the way the stream starsare released from their progenitor. This will be even moreadvantageous in cases where substructure from radial oscil-lations is detectable in the stream.Among those methods that actually model tidal streamsfor comparison to data, N -body simulations are the mostaccurate way to treat the problem, but also the most expen-sive, as well as in some ways the most difficult to interpret.Therefore their application has been limited to date (Howleyet al. 2008; Fardal et al. 2013). Simpler methods generate asingle track to be fitted to the stream. Many authors havesimply assumed that the stream follows the orbit of the pro-genitor, which is a useful approximation but not a generallycorrect one. Other methods for generating a stream trackinclude the streakline techniques of Varghese et al. (2011)and K¨upper et al. (2012), which release particles at constantintervals. These methods are most useful when the progen-itor is slowly and continuously disrupting. Other, more ap-proximate methods for estimating the stream track which c (cid:13) a r X i v : . [ a s t r o - ph . GA ] O c t M. A. Fardal et al. are useful in different circumstances have been used as well(Johnston 1998; Fardal et al. 2006).Action-angle methods, as used by numerous authors(Helmi & White 1999; Eyre & Binney 2011; Sanders & Bin-ney 2013a,b; Bovy 2014; Sanders 2014) result in the clearestconceptual picture of the dynamical structure of the stream.However, these methods are equivalent in their function todoing orbit integrations in the host potential, and there-fore still require accurate release conditions for the particles.Some recipes have been proposed within this context (Eyre& Binney 2011; Bovy 2014), but they have not yet beenshown to be precise or generally applicable. A separate is-sue is that computation of action-angle coordinates requiresseparable potentials for maximum efficiency, and if that isnot possible actions and angles may not exist globally forthe given potential.Recently, Gibbons et al. (2014) and Bonaca et al. (2014)have both proposed modeling the streams using sprays ofparticles released from the satellite and integrated withinthe host potential. Gibbons et al. (2014) assert the poten-tial of the satellite must be included to make this feasible.Using somewhat different release conditions, Bonaca et al.(2014) find it unnecessary to include the satellite potential.However, the testing presented for their method is quite lim-ited and one may question its general applicability. Boththese prescriptions also assume a constant progenitor mass,despite the fact that stream creation implies an ongoing lossof mass.In this paper, we test a recipe for stream particle ejec-tion against a sample of N -body simulations. We judge ourrecipe for release conditions by whether it reproduces distri-butions of actions and frequencies computed from the sim-ulations, since these are the quantities most relevant to theobserved structure of the stream. To make this comparisonpractical, we conduct all the test simulations in an isochronepotential, since the actions and frequencies in this poten-tial are analytic functions of the positions and velocities. Tokeep the testing simple, our progenitor satellites are simpleone-component models such as might be expected for globu-lar clusters. The recipe can be generalized to more complexcases with minor changes.In Section 2, we review the basic properties of tidalstreams, and define the parameters controlling the waystream stars are released. We describe the sample of sim-ulations in Section 3. In Section 4, we constrain the massloss recipe which determines how many stream stars are cre-ated and at what times, as well as the release parameters,which determine the orbits of the stars once created. Sec-tion 5 illustrates the performance of the generated streams.In Section 6, we discuss the wider applicability of our resultsand compare to some other recent stream modeling, and wesummarize our results in Section 7. Very recently, N. Amorisco posted a preprint covering some ofthe issues discussed here (astro-ph/1410.0360).
Although descriptions in other coordinate bases are oftenpossible, the most universal way of describing a tidal streamis as a collection of tracer particles thrown off a satellite asit orbits its host potential, with those particles integratedin position-velocity space (phase space) under the combinedinfluence of the host and the satellite. In the case of a circu-lar orbit within a spherical potential, a necessary conditionfor particles to escape from a satellite at radius r is thatthey exceed the tidal radius (or Jacobi radius), which canbe written as r tidal = (cid:18) m sat f M ( < r ) (cid:19) / r (1)where M is the mass of the host in which the progenitororbits, and f ≡ − Ω − d Φ /dr = 3 − d ln( M ) /d ln( r ) (2)and Ω = V c ( r ) /r = ( GM ( < r ) /r ) / is the circular angu-lar speed (Binney & Tremaine 2008). Note that f = 3 fora Kepler potential and f = 2 for a logarithmic halo. Weassume m sat ( < r tidal ) ≈ m sat , which is commonly satisfiedat least after tides have had a chance to strip the outer-most particles. Even in this case the escape process is notsimple, and much attention has been devoted to the orbitsof stars in the classic point-mass Hill problem and its gen-eralizations (Heggie 2001a). Stars escape into leading andtrailing streams from near the Lagrange points located at ± r tidal from the satellite along the radial vector to the hostcenter.For an eccentric orbit, there is no known rigorous def-inition of a tidal radius. but the concept is still useful inestimating the amount of mass loss. We continue to definethe tidal radius by equation 1, where Ω = V c ( r ) /r is thecircular angular speed. Thus the generalized tidal radius de-pends only on the potential, satellite mass, and radius. Notethat some authors cited here instead set Ω in equation 2 tobe the actual angular velocity v t /r of the satellite, where v t is the tangential velocity.If the stripping is not too violent, the stars again escapeinto well-separated leading and trailing streams from nearthe Lagrange points. A greater eccentricity of the orbit leadsto a disruption rate more strongly peaked near pericenter.If the tidal forces are strong enough, the entire satellite maydisrupt at once, giving a continuum of ejected particles fromleading to trailing. The simplest description of tidal streams occurs when thestream stars can be described by action-angle variables(Helmi & White 1999; Eyre & Binney 2011; Bovy 2014).In these coordinates the stars behave as free particles, atleast if we ignore the continued effect of the satellite as inthe previous subsection. There are certainly cases where anaction-angle description is not possible. For example, thestars may become unbound from the host potential, andhence cease to execute oscillatory behavior. When the po-tential is irregular, as is true for almost any realistic galaxy c (cid:13)000
Although descriptions in other coordinate bases are oftenpossible, the most universal way of describing a tidal streamis as a collection of tracer particles thrown off a satellite asit orbits its host potential, with those particles integratedin position-velocity space (phase space) under the combinedinfluence of the host and the satellite. In the case of a circu-lar orbit within a spherical potential, a necessary conditionfor particles to escape from a satellite at radius r is thatthey exceed the tidal radius (or Jacobi radius), which canbe written as r tidal = (cid:18) m sat f M ( < r ) (cid:19) / r (1)where M is the mass of the host in which the progenitororbits, and f ≡ − Ω − d Φ /dr = 3 − d ln( M ) /d ln( r ) (2)and Ω = V c ( r ) /r = ( GM ( < r ) /r ) / is the circular angu-lar speed (Binney & Tremaine 2008). Note that f = 3 fora Kepler potential and f = 2 for a logarithmic halo. Weassume m sat ( < r tidal ) ≈ m sat , which is commonly satisfiedat least after tides have had a chance to strip the outer-most particles. Even in this case the escape process is notsimple, and much attention has been devoted to the orbitsof stars in the classic point-mass Hill problem and its gen-eralizations (Heggie 2001a). Stars escape into leading andtrailing streams from near the Lagrange points located at ± r tidal from the satellite along the radial vector to the hostcenter.For an eccentric orbit, there is no known rigorous def-inition of a tidal radius. but the concept is still useful inestimating the amount of mass loss. We continue to definethe tidal radius by equation 1, where Ω = V c ( r ) /r is thecircular angular speed. Thus the generalized tidal radius de-pends only on the potential, satellite mass, and radius. Notethat some authors cited here instead set Ω in equation 2 tobe the actual angular velocity v t /r of the satellite, where v t is the tangential velocity.If the stripping is not too violent, the stars again escapeinto well-separated leading and trailing streams from nearthe Lagrange points. A greater eccentricity of the orbit leadsto a disruption rate more strongly peaked near pericenter.If the tidal forces are strong enough, the entire satellite maydisrupt at once, giving a continuum of ejected particles fromleading to trailing. The simplest description of tidal streams occurs when thestream stars can be described by action-angle variables(Helmi & White 1999; Eyre & Binney 2011; Bovy 2014).In these coordinates the stars behave as free particles, atleast if we ignore the continued effect of the satellite as inthe previous subsection. There are certainly cases where anaction-angle description is not possible. For example, thestars may become unbound from the host potential, andhence cease to execute oscillatory behavior. When the po-tential is irregular, as is true for almost any realistic galaxy c (cid:13)000 , 1–20 eneration of mock tidal streams potential, actions do not truly exist. However, the action-angle description continues to be useful, if only as an ap-proximation. For example, in a spherical halo whose poten-tial is flattened by the disk, an action-angle computation islikely to be effective as long as the orbital plane precessionfrequency remains much smaller than either of the radial orazimuthal frequencies.Here we summarize the action-angle approach, mostlyadopting the notation of Bovy (2014). We take the cur-rent time at which the stream is observed to be t = 0.Stars are released from the progenitor at positive lookbacktimes t s , (times t = − t s ). The Hamiltonian is expressedin terms of the actions and is independent of the angles.Therefore the actions are constant with time, and the gra-dient of the Hamiltonian in action space gives the timederivative of the angles, ˙ θ = ∂H/∂ J = Ω ( t s ) = constant,implying θ = θ ( t s ) + Ω t s . With reference to the pro-genitor’s quantities θ ( p ) , Ω ( p ) , the offsets are measured as∆ θ = ∆ θ ( − t s ) + ∆ Ω t s , ∆ Ω = constant. Normally, we canassume that ∆ Ω is small (comparable to r t /r , so that alinear expansion of the frequencies is applicable. This willbe adequate unless the progenitor mass exceeds say 10 − ofits host mass (Bovy 2014), in which case nonlinear termscan become important. Normally, ∆ θ ( t s ) is small as well,so that the properties of ∆ θ are dominated by the secondterm. However, in rare cases, the orbit may be close enoughto circular so that the emitted stars are “trapped” near peri-center and apocenter, rather than following the phase of theprogenitor. We will return to this case later.The orbital behavior of the stars must be coupled witha description of the rate at which particles are released intothe stream, and their initial action and angle offsets. Someauthors focus on models with a continuous release of par-ticles from an initial time to the present (we call this the“leak case”). If the mean frequency of particles released intothe stream is Ω m , then this forms a stream track governedby ∆ θ ≈ ∆ Ω m t s where t s takes on a range of values. Someauthors instead focus on streams formed instantaneously ata single disruption time t d (“burst case”), in which case thetrack is instead described by ∆ θ ≈ t d ∆ Ω where here ∆ Ω takes on a range of values. The formalism of Bovy (2014)has a smooth disruption starting at an finite lookback time t d , so the stream track derived there smoothly interpolatesbetween these two limits at a scale ∆ θ ∼ ∆ Ω m t d , resem-bling the leak case at smaller angles and the burst case atlarger ones.Unless the orbit is nearly circular, we expect the starsto be released predominantly but not entirely near pericen-ter. If the progenitor manages to survive its first pericenterintact, the release model should actually resemble a combi-nation of the leak case with a series of overlaid bursts. Theejection from some satellites might have finite starting andending times, though in this case we expect the scale of thefrequency offsets to change throughout the lifetime of theprogenitor.Although the distribution in action-action space showsa complex bow-tie shape (Eyre & Binney 2011), the distri-bution in frequency-frequency space is nearly diagonal (seeBovy 2014). The distance along the diagonal is correlatedgenerally speaking with the magnitude of the action offset(or more specifically, nearly with the particle energy; John-ston 1998). This leads to two qualitatively different types of “tilts” between the stream and the progenitor orbit. First,for particles released in a single burst, the action and fre-quency offset are strongly correlated. Particles with high en-ergy lag behind in the orbit, while particles with low energyspeed ahead. This creates a tilt between the stream and theorbit.The second kind of tilt results from the fact that thenarrow diagonal distribution of particles in frequency spaceis in most cases not pointed exactly along the frequency vec-tor of the satellite (Eyre & Binney 2011), though in generalit is at least close. One factor in this misalignment angleis that for typical galactic potentials the ratio of radial toazimuthal frequencies tends to decrease for stars executinglarger-radius, higher-energy orbits. Therefore, regardless ofwhether particles are released in a burst or leak slowly, theangle vectors tends to extend along a line misaligned fromthe frequency vector.The action-angle and angle-angle tilts are generallyboth in operation, so that it is not always simple to guesswhere the stream will lie in real space relative to the orbitwithout detailed calculations. However, in young streams(less than a few orbits), we can usually expect the action-angle tilt to dominate the offset between the stream and theorbit. The opposite is true for old streams, say those thathave executed tens of orbits. Naturally, an ejection recipethat correctly reproduces the distribution in release timeand action space should correctly reproduce the effect ofboth of these tilts. Orbits of galactic satellites and globular clusters are typi-cally highly eccentric. Ratios of apocenter r a to pericenter r p of about 4 are typical, as expected for mildly radial tracerpopulations in typical halo potentials (van den Bosch et al.1999). Since tidal forces are a major influence on the massloss rate, and these forces depend on the radius to a highpower ( r − for Kepler or r − for a logarithmic halo), wecan expect the mass loss to proceed as a series of bursts, atleast for the progenitors that do not fall apart in the firstencounter with their host. Thus most tidal streams can beexpected to fall somewhere between the extremes of “burst”and “leak” behavior. This is one possible source of substruc-ture within the stream.K¨upper et al. (2008) and Just et al. (2009) drew at-tention to the substructure in tidal streams. These modelsused progenitors on circular orbits, but subsequent papers(K¨upper et al. 2010, 2012) involved eccentric orbits. A pointpossibly obscured in these papers is that the origin of sub-structure in these two cases is fundamentally different. In thecircular case, the substructure arises from clumps of stars atparticular values of their current radial angle, i.e. the peri-centric locations of the stream stars. This is only possible be-cause the circular case is in the “trapped” regime mentionedabove: the phase of the emitted stars in the trailing (leading)tail is always near pericenter (apocenter), rather than con-tinuously increasing following the phase of the progenitor. Inthe more general case, the substructure is characterized byparticular values of the radial angle at emission , not its cur-rent value. This behavior is clearest in action-angle space,where the stars move in their tori with constant angular ve-locity. In real space, the streaky features (as do all parts of c (cid:13) , 1–20 M. A. Fardal et al. the stream) have their density enhanced at apocenter, dueto the slower motion of the stream particles there. The den-sity seen in the phase-space stream is a combination of thesetwo effects.The dividing line between the trapped and untrappedregimes occurs when the stars emitted at apocenter areplaced on circular orbits. The conditions for this can befound using the epicyclic approximation. We write the az-imuthal frequency as Ω and the radial frequency as κ . Letus suppose that stars are ejected at a position offset by k r r tidal from the satellite, with a tangential velocity offsetby k vφ Ω r tidal from the satellite. To place the ejected starson circular orbits, the radial epicycle amplitude X of thesatellite must take the value X = (2Ω /κ )[ k vφ − (1 + d ln Ω /d ln R ) k r ] r tidal , neglecting second-order terms. For aflat rotation curve this simplifies to X = k vφ r tidal . In otherwords, the epicyclic radius dividing trapped from untrappedregimes is of order r tidal , with the exact value depending onthe values of k r and k vφ (which will be investigated in sec-tion 2.4 below) and the rotation curve. This means thatonly nearly circular orbits (with r a /r p ≈ X ) can bein the trapped regime. For globular clusters with a mass ∼ − that of their host, this demands an epicyclic ra-dius ∼ − that of the orbital radius, or a radial action J r < − L z . For representative dwarf galaxies, this crite-rion rises to J r < − L z . Thus the trapped regime probesa negligibly small region of phase space.In the more common untrapped case which occurs forlarger eccentricity, the clumps seen in simulations spacedfairly evenly in phase along the stream cannot be due to theircurrent radial phase. Instead, they are determined by theradial phase at ejection, and have two simultaneous causes:the variation of the mass loss rate with radial phase, and thecorrelation of action and frequency values with radial phase.The following example illustrates these stream features.Figure 1 (upper panel) shows the particle locations forone of the simulations described in Section 3. This simula-tion is shown at the fourth apocenter passage. We can seestreaky features within the stream, which increasingly over-lap as we go outwards from the satellite. If we convert toaction-angle space (lower panel), the origin of the streaks isclearer. We can see that the streaks correspond to bursts ofparticles released near pericenter, as indicated by the colorcoding. Even though one might expect apocentric clumps ofparticles, due to the long time that the progenitor spendsnear apocenter, these clumps do not appear to exist simplybecause few particles are released around apocenter. Clearly,one must take both the variation in release rate and the dis-persion in the release conditions into account to model thestreaks.In action-angle space, the inclination of the streakshas a consistent sign across both the leading and trailingstreams. The tilt of the streaky features in action-anglespace is ∆Ω = t − s ∆ θ , thus going to zero for large lookbacktimes. This tilt pattern translates fairly well to the real-space view as well. In general streaky features are inclinedsuch that their leading part points inward and their trail-ing part outwards when compared to the mean stream path.Thus streaky features provide an indicator of the stream’sdirection of motion, which may be useful in cases where itis not instantly apparent from other characteristics such as
20 15 10 5 0 5 10 15 20X (kpc)50510152025 Y ( k p c ) R a d i a l f r e q u e n c y Figure 1.
Top: particle positions in the X - Y plane in the orb1.5 m6.0 fat run, described in section 3 and Table 1. Onlyunbound particles are shown, at the fourth apocenter passage ofthe satellite. Color indicates the cosine of the estimated releasephase of the satellite, with red marking pericenter and blue apoc-enter. Bottom: radial frequency versus radial angle, shown at thesame time as the upper panel. kinematic measurements or leading-trailing stream asymme-tries.Counting the streaky features provides a simple in-dication of the number of pericentric passages the pro-genitor has undergone, and thus the timescale over whichthe stream stars have been disrupted. The spacing be-tween the midpoint of the streaky features in angle space is∆Ω m ( r peri ) T radial , constant in time, though the visibility de-creases as time goes on due to the increasing alignment andoverlap of the features. In real space, this spacing will fluc-tuate as the stream moves from apocenter to pericenter andits angular speed varies. The stream tail at | ∆ θ | > | ∆Ω m | t d discussed by Bovy (2014) is in essence the oldest streaky fea-ture in the stream, caused by the burst of particles releasedat the initial pericenter. It would be useful to be able to model the ejection fromthe satellite without performing a full N -body simulation.To model the location of the stream correctly, we need anaccurate recipe for ejecting particles from the satellite. This c (cid:13)000
Top: particle positions in the X - Y plane in the orb1.5 m6.0 fat run, described in section 3 and Table 1. Onlyunbound particles are shown, at the fourth apocenter passage ofthe satellite. Color indicates the cosine of the estimated releasephase of the satellite, with red marking pericenter and blue apoc-enter. Bottom: radial frequency versus radial angle, shown at thesame time as the upper panel. kinematic measurements or leading-trailing stream asymme-tries.Counting the streaky features provides a simple in-dication of the number of pericentric passages the pro-genitor has undergone, and thus the timescale over whichthe stream stars have been disrupted. The spacing be-tween the midpoint of the streaky features in angle space is∆Ω m ( r peri ) T radial , constant in time, though the visibility de-creases as time goes on due to the increasing alignment andoverlap of the features. In real space, this spacing will fluc-tuate as the stream moves from apocenter to pericenter andits angular speed varies. The stream tail at | ∆ θ | > | ∆Ω m | t d discussed by Bovy (2014) is in essence the oldest streaky fea-ture in the stream, caused by the burst of particles releasedat the initial pericenter. It would be useful to be able to model the ejection fromthe satellite without performing a full N -body simulation.To model the location of the stream correctly, we need anaccurate recipe for ejecting particles from the satellite. This c (cid:13)000 , 1–20 eneration of mock tidal streams recipe must supply both the ejection rate from the satelliteand the properties of the ejected particles at the time ofejection. Simple models for the ejection rate could assumea single burst, multiple bursts, or a continuous leak of par-ticles. A more realistic model will assign a distribution ofejection times within a radial cycle.After escaping from the satellite the stars feel its influ-ence diminish rapidly, but the influence of the satellite (andperhaps of the tidal debris) does not stop entirely at thetidal radius for any star (as examined in Choi et al. 2007).This influence is felt especially long for stars that lap theprogenitor in azimuth, or those that have orbital resonanceswith it, and also in cases where the progenitor has a largemass and thus affects even particles at large angular sepa-rations from it. Thus it is not quite correct to examine theconditions of particles crossing the tidal radius and concludethat they determine the final tidal stream structure. How-ever, for small progenitor masses we can still specify thebehavior of stars fairly well in terms of orbits in the hostpotential alone, when ejected from some location near thesatellite with some velocity (not necessarily the actual ve-locity at that point). Let us rewrite these ejection conditionsas follows, using polar coordinates ( r, φ, v r , v t ). r = r sat + k r r tidal (3) φ = φ sat + k φ r tidal /r (4) v r = v r, sat + k vr v r, sat (5) v t = v t, sat + k vφ V c ( r ) r tidal /r (6) z = k z r tidal /r (7) v z = k vz V c ( r sat ) r tidal /r . (8)These forms are similar to some previous models, inparticular the streakline models of Varghese et al. (2011) andK¨upper et al. (2012), The primary difference between thosetwo models is the treatment of the tangential velocity of theejected particles. Varghese et al. (2011) assume a physicalvelocity equal to that of the satellite ( k vφ = 0), while K¨upperet al. (2012) assume an angular velocity equal to that of thesatellite ( k vφ = 1, for the k r = 1 assumed there).In principle the six parameters k r , k φ , k vr , k vφ , k z , k vz could be constants, functions of time, or random samplesfrom distributions that depend on time; we will assume thelatter. Along with other authors we will set k φ = k vr = 0 (9)and test whether this gives a reasonable model for the ejecta.We assume the other parameters are either constants or aredescribed by Gaussian distributions. This is unlikely to beaccurate in the extremes of the distribution, but in this pa-per we are only aiming for an approximation that reproducesthe typical dispersions.The ejection conditions must be coupled with a descrip-tion of the ejection rate, which will be a function of thesatellite’s structure and orbit. With these two ingredients,we can integrate particles in the host potential to predictthe phase-space properties of the stream. To constrain our particle-spray model, we conduct a set of N -body simulations. Our goal is to demonstrate that a rea- Table 1.
Test simulationsName r p (kpc) r a (kpc) r a /r p T rad m sat f t orb1.2 m6.0 .
92 20 .
30 1 . . . orb1.5 m6.0 .
88 22 .
31 1 . . . orb2.0 m6.0 .
40 24 .
77 2 . . . orb3.0 m6.0 .
24 27 .
72 3 . . . orb5.0 m6.0 .
12 30 .
51 5 . . . orb1.2 m6.0 fat .
92 20 .
30 1 . . . orb1.5 m6.0 fat .
88 22 .
31 1 . . . orb2.0 m6.0 fat .
40 24 .
77 2 . . . orb2.0 m6.0 thin .
40 24 .
77 2 . . . orb3.0 m6.0 thin .
24 27 .
72 3 . . . orb5.0 m6.0 thin .
12 30 .
51 5 . . . orb7.0 m6.0 thin .
54 31 .
77 7 . . . sonable model can be derived for some range of satellitesand orbits, but we do not expect that this model will be ap-plicable to all situations. All simulations take place withinan isochrone host potential, so we can calculate actions andangles easily. Our simulations use relatively diffuse satelliteswhere tidal forces are responsible for the disruption. Thisstands in contrast to many models of globular cluster evap-oration, where ejection is often dominated by internal relax-ation processes and stellar mass loss; these must be modeledeither by collisional N -body simulations or Fokker-Planckcodes (Aarseth & Heggie 1998; Chernoff & Weinberg 1990;Murali & Weinberg 1997; Gnedin & Ostriker 1997; Taka-hashi & Portegies Zwart 2000; Lamers et al. 2010; Gieleset al. 2014). Thus the specific fits derived here may not beapplicable to all globular clusters. However, roughly half ofMW clusters (preferentially at small galactocentric radii)have destruction time scales dominated by tidal forces asopposed to relaxation (Heggie 2001b). The observable glob-ular streams may be predominantly those where strong tidalforces induce rapid leakage of stars, so our models are notnecessarily inapplicable to some globulars as well as dwarfgalaxies.Our simulations are of course scale-free, but for com-parison to observed systems we choose a unit system wherethe isochrone potential has mass M = 2 . × M (cid:12) , andscale length b = 3 .
64 kpc. This follows the choices of Eyre& Binney (2011) and provides a rotation velocity similar tothat of the Milky Way for radii near that of the Sun. Allorbits are chosen to have orbital periods 400 Gyr, and the r a /r p span a range 1.2–7.0. We use spherical King (1966)model satellites with W = 3 (i.e., not very concentrated).The simulations are listed in Table 1, and are named accord-ing to the orbit, mass, and scale length.We place the model satellites at the apocenter of theirorbit. We set the initial outer or “tidal” radius of the Kingmodel to be a specified fraction f t of the tidal radius at apocenter : f t ≡ r King /r tidal ( r apo ) . (10)The tidal radius at pericenter can be inferred from Table 1and equation 1, but in most cases it is inside the initialmodel outer radius, inducing significant mass loss at eachpericenter.We run the test simulations until the time reaches 4.3Gyr. We save snapshots every 10 Myr. We track the positionof the satellite as a function of time and convert it to action- c (cid:13) , 1–20 M. A. Fardal et al. orb1.2_m6.0 orb5.0_m6.0 orb2.0_m6.0_thin orb1.5_m6.0 orb1.2_m6.0_fat orb3.0_m6.0_thin orb2.0_m6.0 orb1.5_m6.0_fat orb5.0_m6.0_thin orb3.0_m6.0 orb2.0_m6.0_fat orb7.0_m6.0_thinEjection time (Gyr) S a t e lli t e m a ss ( M fl ) Figure 2.
Satellite mass as a function of time. Black: mass within tidal radius measured from simulations. Wiggles illustrating periodicmass loss enhancements are clearly visible. Red: predicted mass at successive apocenters, obtained by truncating actual simulationstructure at the previous apocenter near the pericentric tidal radius, as discussed in the text. Green: mass predicted from the full massloss model, which uses only the initial satellite orbit, mass and structure as inputs. Small symbols denote mass at apocenter. angle variables. In some cases the orbital frequencies seemto deviate slightly from the expectations based on the initialorbit, most likely due to interaction of the satellite with itsown tidal debris (as expected from Choi et al. 2007). Thus wecompute corrected frequencies from the measured satellitepositions, and use the new frequencies to initialize a finelyspaced lookup table for each satellite’s position and velocityas a function of simulation time.For each particle outside the tidal radius at the finalsnapshot, we estimate the ejection time from its action-anglecoordinates. Using this to estimate a useful starting time, wecompute the orbit of the particle backwards in time betweensnapshots in a two-body time-dependent potential, using thelookup table to compute the satellite position. We use thisto infer a more accurate ejection time from the satellite, asdefined by the time at which the particle crosses the tidalradius. We then save the position and velocity of the par-ticle and satellite. In some cases this procedure fails. Forsome cases, the radial and azimuthal angles do not predictthe same ejection time. Many of these are near the satel-lite and some have probably interacted multiple times withthe satellite. Some other failures occur when the satellite isfalling apart rapidly and the point mass approximation fails.When the refined approximation fails, we find the ejectiontime using a linear interpolation of the radial distance fromthe satellite at snapshots.
To form a complete model of the tidal stream, we mustmodel the rate of mass loss in addition to the conditions ofthe ejected particles. We will do this by comparison to thesimulations in section 3. The first step in constructing themass model is to specify the satellite mass at subsequentapocenters. Earlier work on satellites of various types hassuggested that the evolution of the satellite mass and den-sity profile follow quite predictable tracks, given only theinitial orbit and structure (Hayashi et al. 2003; Pe˜narrubiaet al. 2008). In our simulation set, we find that the satellitemass at the next apocenter can be reproduced quite well inmost cases simply by truncating the apocentric profile at 0.9times the tidal radius at pericenter. Figure 2 illustrates theresults of this procedure. The red curve represents the satel-lite mass predicted at each apocenter, based on truncatingthe satellite structure of the previous apocenter. This agreeswell with the simulation results shown by the black curve.However, it is not sufficient for our model to predict themass using the actual, evolving structure of the simulated N -body satellite, as we want to dispense with N -body sim-ulations altogether. We find that after the first pericentricpasssage, the radial profiles of satellites that do not disruptaltogether are reasonably well described by Einasto profiles, c (cid:13)000
To form a complete model of the tidal stream, we mustmodel the rate of mass loss in addition to the conditions ofthe ejected particles. We will do this by comparison to thesimulations in section 3. The first step in constructing themass model is to specify the satellite mass at subsequentapocenters. Earlier work on satellites of various types hassuggested that the evolution of the satellite mass and den-sity profile follow quite predictable tracks, given only theinitial orbit and structure (Hayashi et al. 2003; Pe˜narrubiaet al. 2008). In our simulation set, we find that the satellitemass at the next apocenter can be reproduced quite well inmost cases simply by truncating the apocentric profile at 0.9times the tidal radius at pericenter. Figure 2 illustrates theresults of this procedure. The red curve represents the satel-lite mass predicted at each apocenter, based on truncatingthe satellite structure of the previous apocenter. This agreeswell with the simulation results shown by the black curve.However, it is not sufficient for our model to predict themass using the actual, evolving structure of the simulated N -body satellite, as we want to dispense with N -body sim-ulations altogether. We find that after the first pericentricpasssage, the radial profiles of satellites that do not disruptaltogether are reasonably well described by Einasto profiles, c (cid:13)000 , 1–20 eneration of mock tidal streams orb1.2_m6.0 orb5.0_m6.0 orb2.0_m6.0_thin orb1.5_m6.0 orb1.2_m6.0_fat orb3.0_m6.0_thin orb2.0_m6.0 orb1.5_m6.0_fat orb5.0_m6.0_thin -2 -1 orb3.0_m6.0 -2 -1 orb2.0_m6.0_fat -2 -1 orb7.0_m6.0_thinRadius (kpc) D e n s i t y ( M fl k p c − ) Figure 3.
Radial density profiles, as measured at the third apocenter. Black: profile measured from simulation, discarding any formallyunbound particles. Red: profile generated from mass loss model using current mass, initial mass, and initial density profile. Blue: tidalradius at apocenter. as shown in Figure 3. In this figure, we only include thoseparticles that are formally bound to the satellite, as definedby the relative velocity and potential from the satellite andits debris (and ignoring the tidal force). We write the Einastoprofile in the form ρ ( r ) = (2 µ ) µ m sat πr sc µe µ Γ(3 µ ) exp (cid:26) − µ (cid:20)(cid:16) rr sc (cid:17) /µ − (cid:21)(cid:27) (11)so the index µ resembles the n in the related Sersic profile.Here r sc is the radius at which the logarithmic slope of thedensity equals the isothermal slope −
2. For this profile, thecumulative mass is given by M ( < r ) = m sat γ (cid:2) µ, µ ( r/r sc ) /µ (cid:3) (12)where γ is the normalized incomplete gamma function.Furthermore, we find that the radial scale evolves littleover the course of the simulation, while the density scaledrops significantly. Therefore we simply assume the radialscale is fixed at r sc = 0 . r outer , where r outer is the initialouter or “tidal” radius of the King model. The density scalethen simply scales with the initial mass. The optimal µ variesslightly between runs and snapshots, but for simplicity wekeep it fixed at µ = 0 .
9. In other words, the density profileis fairly close to exponential.The resulting fits are shown in Figure 3. At large radiusthe tidal debris shows some complicated structure whichcannot be modeled by a simple profile, but within a radiuscontaining most of the mass the fit is fairly good. Once we have a sequence of apocenter masses, this determines theamount of mass lost over each radial cycle, and this predic-tion is shown in Figure 2.As our simulations take place in a spherical potential,we include only the effects of “bulge shocking” but neglectthose of “disk shocking”. The performance of our recipeshould be checked particularly in cases where satellites passdirectly through the disk of the host, as these may requirean enhanced degree of mass loss which is dependent on thefull details of the orbit rather than just the oscillations inradius.To determine the variation of mass ejection rate withinthis cycle, we first define the “acceleration gradient” as g a ≡ Ω − d Φ( r ) /dr , (13)where Ω = V c ( r ) /r is again the circular frequency. Thisquantity represents the second derivative with radius of theeffective potential, for a circular orbit, and thus measuresthe strength of the tidal force. (Note g a = Ω f ). We thencompute the ratio of this quantity at apocenter and pericen-ter, R acc ≡ g a ( r peri ) g a ( r apo ) (14)We use R acc as our measure for the variation of the tidalforces over a radial oscillation.We use the following analytic form to define the ejectionrate as a function of radial phase θ r : c (cid:13) , 1–20 M. A. Fardal et al. orb1.2_m6.0 orb5.0_m6.0 orb2.0_m6.0_thin orb1.5_m6.0 orb1.2_m6.0_fat orb3.0_m6.0_thin orb2.0_m6.0 orb1.5_m6.0_fat orb5.0_m6.0_thin orb3.0_m6.0 orb2.0_m6.0_fat orb7.0_m6.0_thinRadial phase at ejection θ r E j e c t i o n p h a s e d i s t r i b u t i o n ( n o r m a li z e d ) Figure 4.
Distribution of mass loss in ejection radial phase. Blue histograms: distribution measured from simulations, as measuredfrom crossing at one tidal radius. Red: model specified by orbit and initial size, as in text. Clearly the model captures the overall peakybehavior while ignoring some of the finer structure. dMdθ r = A (cid:18) r ej − (cid:20) θ r − θ mid )2 (cid:21) α (cid:19) (15)With this form, the ratio of ejection at peak to trough is r ej . The normalization and cumulative distribution of thisform can be calculated using beta functions (complete andincomplete respectively).We then set the parameters in this expression using r ej = exp(1 . R / acc ) (16) α = R . acc (17) θ mid = − . . f t R acc f t R acc (18)These forms were determined by comparison with the dis-tribution of ejection phase in the simulations (Figure 4).These parameter expressions are sufficient to approx-imate the variation of emission rate with satellite orbitalphase, to the degree required to specify reasonable models.They are not well justified physically so should be checkedin cases that are far from Note that the rate of particle emis-sion peaks substantially after pericenter in some cases.With these specifications, we can now model the parti-cle emission from the simulated satellite. We first determinethe number of particles to be ejected between subsequentapocenters, then randomly sample from equation 15 usingthe rejection method, which generates a sequence of ejectedparticles. The mass in the satellite is determined by the an-alytic formula for the cumulative ejected mass. The last step is to specify the orbital release conditions of theejected particles. We determine this by comparison with thefinal actions and frequencies of the ejected particles in thesimulations as a function of radial phase at ejection. Sincewe are using spherical isochrone potential, the azimuthaland vertical frequencies are degenerate, so we focus on thebehavior in 4-dimensional action-angle space.We specify the release using equations 3–6, setting k φ = 0, k vr = 0, and choosing k r , k vφ to be constants. Webegan by examining the predictions of the Varghese et al.(2011) and K¨upper et al. (2012) initial conditions discussedpreviously. The mean action and frequency offsets generatedby the two initial conditions are fairly similar when averagedover a radial cycle, but the variation of the frequencies withphase is much stronger using the K¨upper initial conditions.We optimized the choice of constants essentially by eye. Themean values for the k constants are¯ k r = 2 . , (19)¯ k vφ = 0 . . (20)This is closer to the Varghese conditions than the K¨uppercondition, the choice of which tends to make the frequenciesvary too much with orbital phase. The optimal parametersappear to vary slightly between simulations, but not stronglyenough for us to justify a more complex choice.However, Figures 5 and 6 show that the dispersion in c (cid:13)000
Distribution of mass loss in ejection radial phase. Blue histograms: distribution measured from simulations, as measuredfrom crossing at one tidal radius. Red: model specified by orbit and initial size, as in text. Clearly the model captures the overall peakybehavior while ignoring some of the finer structure. dMdθ r = A (cid:18) r ej − (cid:20) θ r − θ mid )2 (cid:21) α (cid:19) (15)With this form, the ratio of ejection at peak to trough is r ej . The normalization and cumulative distribution of thisform can be calculated using beta functions (complete andincomplete respectively).We then set the parameters in this expression using r ej = exp(1 . R / acc ) (16) α = R . acc (17) θ mid = − . . f t R acc f t R acc (18)These forms were determined by comparison with the dis-tribution of ejection phase in the simulations (Figure 4).These parameter expressions are sufficient to approx-imate the variation of emission rate with satellite orbitalphase, to the degree required to specify reasonable models.They are not well justified physically so should be checkedin cases that are far from Note that the rate of particle emis-sion peaks substantially after pericenter in some cases.With these specifications, we can now model the parti-cle emission from the simulated satellite. We first determinethe number of particles to be ejected between subsequentapocenters, then randomly sample from equation 15 usingthe rejection method, which generates a sequence of ejectedparticles. The mass in the satellite is determined by the an-alytic formula for the cumulative ejected mass. The last step is to specify the orbital release conditions of theejected particles. We determine this by comparison with thefinal actions and frequencies of the ejected particles in thesimulations as a function of radial phase at ejection. Sincewe are using spherical isochrone potential, the azimuthaland vertical frequencies are degenerate, so we focus on thebehavior in 4-dimensional action-angle space.We specify the release using equations 3–6, setting k φ = 0, k vr = 0, and choosing k r , k vφ to be constants. Webegan by examining the predictions of the Varghese et al.(2011) and K¨upper et al. (2012) initial conditions discussedpreviously. The mean action and frequency offsets generatedby the two initial conditions are fairly similar when averagedover a radial cycle, but the variation of the frequencies withphase is much stronger using the K¨upper initial conditions.We optimized the choice of constants essentially by eye. Themean values for the k constants are¯ k r = 2 . , (19)¯ k vφ = 0 . . (20)This is closer to the Varghese conditions than the K¨uppercondition, the choice of which tends to make the frequenciesvary too much with orbital phase. The optimal parametersappear to vary slightly between simulations, but not stronglyenough for us to justify a more complex choice.However, Figures 5 and 6 show that the dispersion in c (cid:13)000 , 1–20 eneration of mock tidal streams N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z orb1.2_m6.0orb1.5_m6.0orb2.0_m6.0orb3.0_m6.0 Ejection phase Figure 5.
Actions of stream stars. First two columns show J r , second two show L z . Both are given in units of kpc km s − . In each pairof columns, simulation results are on the left and results from the mass loss model on the right. N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z orb5.0_m6.0orb1.2_m6.0_fatorb1.5_m6.0_fatorb2.0_m6.0_fat Ejection phase Figure 5 – continued c (cid:13) , 1–20 M. A. Fardal et al. N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ J r M o c k s t r e a m ∆ J r N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z N - b o d y ∆ L z M o c k s t r e a m ∆ L z orb2.0_m6.0_thinorb3.0_m6.0_thinorb5.0_m6.0_thinorb7.0_m6.0_thin Ejection phase Figure 5 – continued N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ orb1.2_m6.0orb1.5_m6.0orb2.0_m6.0orb3.0_m6.0 Ejection phase Figure 6.
Orbital frequencies of stream stars. First two columns show Ω r , second two show Ω φ . These are given in units of radians /Gyr. In each pair of columns, simulation results are on the left and results from the mass loss model on the right.c (cid:13)000
Orbital frequencies of stream stars. First two columns show Ω r , second two show Ω φ . These are given in units of radians /Gyr. In each pair of columns, simulation results are on the left and results from the mass loss model on the right.c (cid:13)000 , 1–20 eneration of mock tidal streams N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ orb5.0_m6.0orb1.2_m6.0_fatorb1.5_m6.0_fatorb2.0_m6.0_fat Ejection phase Figure 6 – continued N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω r M o c k s t r e a m ∆ Ω r N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ N - b o d y ∆ Ω φ M o c k s t r e a m ∆ Ω φ orb2.0_m6.0_thinorb3.0_m6.0_thinorb5.0_m6.0_thinorb7.0_m6.0_thin Ejection phase Figure 6 – continued c (cid:13) , 1–20 M. A. Fardal et al. L r a d ( k p c k m s − ) ( s i m ) L r a d ( k p c k m s − ) ( m o c k ) L t a n g ( k p c k m s − ) ( s i m ) L t a n g ( k p c k m s − ) ( m o c k ) orb3.0_m6.0 Figure 7.
Angular momentum components in orbital plane, usedas indicators of off-plane motions of the particles, shown as afunction of ejected phase. First two columns show the compo-nent along the radial vector to the satellite, while the second twoshow the component along the tangential vector. In each pair ofcolumns, simulation results are on the left and results from themass loss model on the right. release conditions is not constant, but decreases significantlyfor simulations with little mass loss (those on near-circularorbits). We address this by choosing dispersions σ ( k r ) and σ ( k vφ ) of the form σ ( k r ) = σ ( k vt ) = min(0 . f t R / acc , .
4) (21)The initial conditions for the particle-spray models in theseand subsequent figures include these dispersions.We also introduce a random nonzero offset in the off-plane position and velocity. For spherical simulations, this isnecessary if we are to prevent the model stream from beingentirely flat, in contrast to the simulation results. In a non-spherical potential, the vertical extent also thickens due todifferences between Ω z and Ω φ . We choose σ ( k z ) and σ ( k vz )to reproduce the dispersion in the off-plane components ofangular momentum along the radial and tangential vectorsto the satellite at the time of ejection. In our model givenby equations 7–8, these two components are determined by k z and k vz respectively. The choices¯ k z = ¯ k vz = 0 (by symmetry)) (22) σ ( k z ) = 0 . σ ( k vz ) = 0 . k r , k vφ k z , and k vz from Gaussian distributions withthe specified mean and dispersion. We can then determinethe positions and velocities using equations 3–8 and 19–20and then determine the actions (Figure 5) and frequencies(Figure 6) for each of our ejected particles. These figureswere used to guide our choice of constants. The action plotsof the simulations show a significant amount of interestingstructure in some cases, in contrast to the particle-spraymodel which is fairly symmetrical about the origin and uni-modal at any ejection phase. However, the structure in the frequency plots appears overall less complicated and moresimilar to the particle-spray model, though in some cases bi-modal structures or sprays of particles are still evident. Thefrequencies influence the stream properties more than theactions, so this is an encouraging sign for the particle-spraymodel.The correlation between the angular momentum andradial action of the stream particles can be seen in the secondand fourth rows of Figure 8. Here the first column is thesimulation results, the second is the result of our recipe,while the third and fourth columns will be discussed later.It can be seen that our recipe reproduces the characteristic“bow-tie” pattern noted for this plot in this previous work(Eyre & Binney 2011; Bovy 2014).Once we have chosen ejection times and phase-spacecoordinates for all the mock stream particles, we can evolvethem through time either by direct orbital integration or byusing the action-angle formalism. We use the latter methodhere because of the ease of action-angle calculations in theisochrone potential. We propagate the angles to the finaltimestep using the computed frequencies of each particle.Finally, we can determine the position and velocity of eachparticle by inverting the real-space to action-angle transfor-mation. The results are shown in Figure 9. The agreementis quite impressive overall, showing deviation from the or-bit by the correct amount and a fairly good agreement onthe length, width, and surface brightness of the stream. Themain deviations from the simulation can be found very closeto the satellite and in the extreme tails of the distribution,but a careful inspection is necessary to find these differences. In Figure 8, we compare the results of our recipe to two otherrecipes that have appeared recently, which we implementedto the best of our understanding. The results of the simula-tions and our recipe appear in the first and second columnsrespectively. The Gibbons et al. (2014) recipe, which theycall the “Lagrange Cloud stripping” method, appears in thethird column, though in this case we are not including anadditional force from the satellite as they found necessary,but simply integrating particles in the host potential. Also,here we take the velocity dispersion at particle release to bebe σ v = V c ( r ) r t /r , which we found to be obeyed reasonablyaccurately by measuring the actual 1-d velocity dispersionof bound particles over our ensemble of simulations. In theGibbons approach this dispersion is a free parameter used infitting the results of simulations or observations. The recipeof Bonaca et al. (2014), which they call the “Fast-forward”method, appears in the fourth column. Both these recipesrequire an estimate for the velocity dispersion, but this istaken to be a free parameter used in fitting in the Gibbonsapproach.For the two simulation cases shown, the Gibbons reciperesults in frequency offsets that are generally too small. Thisagrees with their finding that the streams are too shortwithout an additional force from the satellite. In contrast,the scale of the offset in the Bonaca recipe is fairly good. c (cid:13)000
4) (21)The initial conditions for the particle-spray models in theseand subsequent figures include these dispersions.We also introduce a random nonzero offset in the off-plane position and velocity. For spherical simulations, this isnecessary if we are to prevent the model stream from beingentirely flat, in contrast to the simulation results. In a non-spherical potential, the vertical extent also thickens due todifferences between Ω z and Ω φ . We choose σ ( k z ) and σ ( k vz )to reproduce the dispersion in the off-plane components ofangular momentum along the radial and tangential vectorsto the satellite at the time of ejection. In our model givenby equations 7–8, these two components are determined by k z and k vz respectively. The choices¯ k z = ¯ k vz = 0 (by symmetry)) (22) σ ( k z ) = 0 . σ ( k vz ) = 0 . k r , k vφ k z , and k vz from Gaussian distributions withthe specified mean and dispersion. We can then determinethe positions and velocities using equations 3–8 and 19–20and then determine the actions (Figure 5) and frequencies(Figure 6) for each of our ejected particles. These figureswere used to guide our choice of constants. The action plotsof the simulations show a significant amount of interestingstructure in some cases, in contrast to the particle-spraymodel which is fairly symmetrical about the origin and uni-modal at any ejection phase. However, the structure in the frequency plots appears overall less complicated and moresimilar to the particle-spray model, though in some cases bi-modal structures or sprays of particles are still evident. Thefrequencies influence the stream properties more than theactions, so this is an encouraging sign for the particle-spraymodel.The correlation between the angular momentum andradial action of the stream particles can be seen in the secondand fourth rows of Figure 8. Here the first column is thesimulation results, the second is the result of our recipe,while the third and fourth columns will be discussed later.It can be seen that our recipe reproduces the characteristic“bow-tie” pattern noted for this plot in this previous work(Eyre & Binney 2011; Bovy 2014).Once we have chosen ejection times and phase-spacecoordinates for all the mock stream particles, we can evolvethem through time either by direct orbital integration or byusing the action-angle formalism. We use the latter methodhere because of the ease of action-angle calculations in theisochrone potential. We propagate the angles to the finaltimestep using the computed frequencies of each particle.Finally, we can determine the position and velocity of eachparticle by inverting the real-space to action-angle transfor-mation. The results are shown in Figure 9. The agreementis quite impressive overall, showing deviation from the or-bit by the correct amount and a fairly good agreement onthe length, width, and surface brightness of the stream. Themain deviations from the simulation can be found very closeto the satellite and in the extreme tails of the distribution,but a careful inspection is necessary to find these differences. In Figure 8, we compare the results of our recipe to two otherrecipes that have appeared recently, which we implementedto the best of our understanding. The results of the simula-tions and our recipe appear in the first and second columnsrespectively. The Gibbons et al. (2014) recipe, which theycall the “Lagrange Cloud stripping” method, appears in thethird column, though in this case we are not including anadditional force from the satellite as they found necessary,but simply integrating particles in the host potential. Also,here we take the velocity dispersion at particle release to bebe σ v = V c ( r ) r t /r , which we found to be obeyed reasonablyaccurately by measuring the actual 1-d velocity dispersionof bound particles over our ensemble of simulations. In theGibbons approach this dispersion is a free parameter used infitting the results of simulations or observations. The recipeof Bonaca et al. (2014), which they call the “Fast-forward”method, appears in the fourth column. Both these recipesrequire an estimate for the velocity dispersion, but this istaken to be a free parameter used in fitting in the Gibbonsapproach.For the two simulation cases shown, the Gibbons reciperesults in frequency offsets that are generally too small. Thisagrees with their finding that the streams are too shortwithout an additional force from the satellite. In contrast,the scale of the offset in the Bonaca recipe is fairly good. c (cid:13)000 , 1–20 eneration of mock tidal streams Ejection phase2.01.51.00.50.00.51.01.52.0 N - b o d y ∆ Ω r orb1.5_m6.0 Ejection phase2.01.51.00.50.00.51.01.52.0 M o c k s t r e a m ∆ Ω r Ejection phase2.01.51.00.50.00.51.01.52.0 G i bb o n s - li k e ∆ Ω r Ejection phase2.01.51.00.50.00.51.01.52.0 B o n a c a - li k e ∆ Ω r L tot N - b o d y J r L tot M o c k s t r e a m J r L tot G i bb o n s - li k e J r L tot B o n a c a - li k e J r Ejection phase42024 N - b o d y ∆ Ω r orb5.0_m6.0 Ejection phase42024 M o c k s t r e a m ∆ Ω r Ejection phase42024 G i bb o n s - li k e ∆ Ω r Ejection phase42024 B o n a c a - li k e ∆ Ω r L tot N - b o d y J r L tot M o c k s t r e a m J r L tot G i bb o n s - li k e J r L tot B o n a c a - li k e J r Figure 8.
Comparison of different algorithms. The first two and second two rows use the orb1.5 m6.0 and orb5.0 m6.0 runs respectively.The first row for each run shows the radial frequency offset versus ejection phase. The second row for each run shows the radial actionversus the angular momentum. The first column shows the simulations, the second our recipe, the third the Gibbons et al. (2014) recipe(minus the force from the satellite), and the fourth the Bonaca et al. (2014) recipe.
However, the contrast between minimum and maximum fre-quency is somewhat too large.The Gibbons streams have frequency distributions thatare too broad, especially for the first run shown ( r a /r p =1 . σ equal to the velocity dispersionof the satellite. In the Gibbons recipe, σ is actually treatedas a free parameter. So in their approach the optimizationshould eventually pick a reasonable dispersion, but at thecost of losing any information obtainable from linking therelease velocity dispersion to the satellite properties.The Bonaca streams in contrast have a large radial ac-tion dispersion, but almost no dispersion in angular momen-tum, owing to the inclusion of dispersion only in the radialvelocity. For certain simulations (see the lowest row in Fig-ure 8), the dispersion is large enough to throw some starsbackwards into the stream emerging from the opposing La-grange point. This results in two extra streams emergingfrom the satellite in these cases. For most simulations, thefrequency dispersion is too small at pericenter, which resultsin streaky features that are too short compared to the simu-lations. In summary, our recipe seems to improve the overallagreement with simulations compared to previous recipes forthe cases examined here, though the differences will be morevisible in some cases (mainly those with particularly small or large eccentricities) than in others. The general similarity,along with the agreement with simulations we demonstratedin section 4, supports the general approach of particle-spraymethods as used here and in the work of Gibbons et al.(2014) and Bonaca et al. (2014).Turning now to the overall ensemble of simulations, themodel streams shown in Figure 9 agree extremely well withthe simulation results. The width, length, and substructurewithin the streams are all quite well reproduced. The mostcommon deficiency is an excessively long tail of debris, whichstems in part from our use of simple Gaussian distributionsof the space and velocity offsets rather than a more sharplytruncated distribution. To extend the test of the method beyond the parameterrange where it has been calibrated, we run an N -body sim-ulation based on the one shown in Figure 2 of Gibbons et al.(2014). This simulation uses a different potential and a muchmore massive satellite than those we have examined so far.The host potential for this run is a spherical Navarro-Frenk-White profile, with a mass of 7 . × M (cid:12) contained withina radius of 185 .
41 kpc, and a scale radius of 9 .
27 kpc yieldinga concentration of 20. This equates to a virial overdensityof 180 with respect to the critical density for a Hubble con-stant of H = 75 km s − . We use a W = 3 King model of6 . × M (cid:12) and outer radius 4 .
75 kpc We set this moving c (cid:13) , 1–20 M. A. Fardal et al. x (kpc)30201001020 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)30201001020 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)30201001020 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c )
30 20 10 0 10 20x (kpc)30201001020 N - b o d y y ( k p c )
30 20 10 0 10 20x (kpc) M o c k s t r e a m y ( k p c ) φ (radians)18012060060120 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)18012060060120 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)18012060060120 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)18012060060120 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r orb1.2_m6.0orb1.5_m6.0orb2.0_m6.0orb3.0_m6.0 Figure 9.
Stream position / velocity. First two columns show a face-on view of the orbit plane, while the second two show the radialvelocity versus the azimuthal angle. In each pair of columns, simulation results are on the left and results from the mass loss model onthe right. Streams are shown at the last apocenter timestep at which the satellite is still intact. x (kpc)20100102030 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)20100102030 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)20100102030 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c )
24 16 8 0 8 16x (kpc)20100102030 N - b o d y y ( k p c )
24 16 8 0 8 16x (kpc) M o c k s t r e a m y ( k p c ) φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r orb5.0_m6.0orb1.2_m6.0_fatorb1.5_m6.0_fatorb2.0_m6.0_fat Figure 9 – continued c (cid:13)000
24 16 8 0 8 16x (kpc) M o c k s t r e a m y ( k p c ) φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)24016080080160 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r orb5.0_m6.0orb1.2_m6.0_fatorb1.5_m6.0_fatorb2.0_m6.0_fat Figure 9 – continued c (cid:13)000 , 1–20 eneration of mock tidal streams x (kpc)45301501530 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)45301501530 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)45301501530 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c )
30 15 0 15 30x (kpc)45301501530 N - b o d y y ( k p c )
30 15 0 15 30x (kpc) M o c k s t r e a m y ( k p c ) φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r orb2.0_m6.0_thinorb3.0_m6.0_thinorb5.0_m6.0_thinorb7.0_m6.0_thin Figure 9 – continued x (kpc)45301501530 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)45301501530 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c ) x (kpc)45301501530 N - b o d y y ( k p c ) x (kpc) M o c k s t r e a m y ( k p c )
30 15 0 15 30x (kpc)45301501530 N - b o d y y ( k p c )
30 15 0 15 30x (kpc) M o c k s t r e a m y ( k p c ) φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r φ (radians)3002001000100200 N - b o d y v r ( k m s − ) φ (radians) M o c k s t r e a m v r orb2.0_m6.0_thinorb3.0_m6.0_thinorb5.0_m6.0_thinorb7.0_m6.0_thin Figure 9 – continued c (cid:13) , 1–20 M. A. Fardal et al.
120 80 40 0 40 80x (kpc)120804004080 N - b o d y y ( k p c )
120 80 40 0 40 80x (kpc)120804004080 M o c k s t r e a m y ( k p c ) φ (radians)765432 N - b o d y E ( k m s − ) φ (radians)765432 M o c k s t r e a m E ( k m s − ) Figure 10. N -body simulation roughly emulating the Sagittarius stream (left column), following Gibbons et al. (2014), and particle-spray model thereof (right column). The top row shows particle positions. The bottom row shows the particle energies versus azimuth;streams from different pericenter passages are particularly well separated in this representation. The progenitor appears at x = 0 kpc,top left vertical streak in The complex structure of the simulated stream is well reproduced by the particle-spray model. in the x - y plane with coordinates x = − .
327 kpc, y =70 .
772 kpc, v x = − .
71 km s − , and v y = − .
59 km s − .The resulting orbit has apocenter 70 . . f t = 0 . .
14 Gyr into the run. At this point thesatellite is just past the fourth pericenter, so three separatestreams from previous pericentric passages are clearly vis-ible in both the leading and trailing streams. The spatialdistribution of the particles is similar to that of Gibbonset al. (2014). The different streams separate particularly wellin the lower row, which shows the particle energies plottedversus azimuth. In the spatial distribution, these streamsseparate most clearly at apocenter of the radial loops. Thissubstructure clearly shows the action-angle tilt which is ex-pected to dominate for a young stream (only 3.5 radial os-cillations), as discussed above.Results from our particle-spray model are shown in theright-hand column. We have performed the required orbitintegrations using the galpy package written by Jo Bovy.Clearly, the complex structure of the stream is well repro-duced by our model. When overlaid, the multiple compo- nents of the leading and trailing streams match the N -bodymodel nearly exactly. Some slight differences are visible, no-tably gaps at two different azimuth values in the leadingstream of the simulation and one in the leading stream inthe energy plot. One of these gaps is also visible in the spatialdistribution, and results in a “C”-shaped feature separatedfrom the rest of the leading stream. These features resultfrom the satellite crossing its own stream, as made inevitableby the perfectly spherical potential used here, and creatinga gap in the manner studied by Carlberg (2013). Since ourmodel does not include the satellite potential, these gapsdo not appear in the right-hand column. The model streamis somewhat longer and in some places narrower than the N -body results. Also, differences are easily apparent nearthe progenitor satellite, both from the progenitor itself (notrepresented in the particle-spray model) and from particlesthat have been stripped but have yet to acquire their full ac-tion offset as they curve away from the progenitor. Overall,however, the agreement is extremely good.Figure 11 shows the performance of the particle sprayrecipe in a compound bulge-disk-halo potential. Two or-bits are shown, one on a regular orbit and one on a mildlychaotic orbit. The satellite in this test is located at a po- c (cid:13)000
14 Gyr into the run. At this point thesatellite is just past the fourth pericenter, so three separatestreams from previous pericentric passages are clearly vis-ible in both the leading and trailing streams. The spatialdistribution of the particles is similar to that of Gibbonset al. (2014). The different streams separate particularly wellin the lower row, which shows the particle energies plottedversus azimuth. In the spatial distribution, these streamsseparate most clearly at apocenter of the radial loops. Thissubstructure clearly shows the action-angle tilt which is ex-pected to dominate for a young stream (only 3.5 radial os-cillations), as discussed above.Results from our particle-spray model are shown in theright-hand column. We have performed the required orbitintegrations using the galpy package written by Jo Bovy.Clearly, the complex structure of the stream is well repro-duced by our model. When overlaid, the multiple compo- nents of the leading and trailing streams match the N -bodymodel nearly exactly. Some slight differences are visible, no-tably gaps at two different azimuth values in the leadingstream of the simulation and one in the leading stream inthe energy plot. One of these gaps is also visible in the spatialdistribution, and results in a “C”-shaped feature separatedfrom the rest of the leading stream. These features resultfrom the satellite crossing its own stream, as made inevitableby the perfectly spherical potential used here, and creatinga gap in the manner studied by Carlberg (2013). Since ourmodel does not include the satellite potential, these gapsdo not appear in the right-hand column. The model streamis somewhat longer and in some places narrower than the N -body results. Also, differences are easily apparent nearthe progenitor satellite, both from the progenitor itself (notrepresented in the particle-spray model) and from particlesthat have been stripped but have yet to acquire their full ac-tion offset as they curve away from the progenitor. Overall,however, the agreement is extremely good.Figure 11 shows the performance of the particle sprayrecipe in a compound bulge-disk-halo potential. Two or-bits are shown, one on a regular orbit and one on a mildlychaotic orbit. The satellite in this test is located at a po- c (cid:13)000 , 1–20 eneration of mock tidal streams
30 25 20 15 10 5x (kpc)252015105 z ( k p c )
30 25 20 15 10 5x (kpc)252015105 z ( k p c )
30 25 20 15 10 5x (kpc)252015105 z ( k p c )
30 25 20 15 10 5x (kpc)252015105 z ( k p c ) Figure 11.
Simulations using a bulge-disk-halo potential, as discussed in the text Left columns: N -body simulations. Right columns:particle-spray models. Upper row: regular orbit with q = 0 .
57. Lower row: apparently chaotic orbit with q = 0 . sition x = − .
59 kpc, y = − .
55 kpc, z = − .
89 kpcand velocity v x = − . − , v y = 24 .
36 km s − , and v z = − .
12 km s − . We initialize the satellite with mass4 × M (cid:12) and tidal factor f t = 0 .
7. The orbit is inte-grated backwards from the given point for t d = 2694 Myr,and then the particle-spray and N -body simulations areevolved to the present day. In the fixed potential used here,the bulge is a Hernquist profile with mass 3 . × M (cid:12) and radius 0 . M (cid:12) , a = 6 . b = 0 .
26 kpc. The halo is anoblate logarithmic halo with φ ( R, z ) = V h ln( R + ( z/q ) + d ), V h = 115 km s − , and d = 12 kpc. The only differencebetween the orbits used in the top and bottom rows is thatthe top uses q = 0 .
57, and the bottom q = 0 . q = 0 .
63, the orbit be-comes regular again, and in tandem the stream becomesnarrow again so that it much more closely resembles the toppanel. (Admittedly, the flattening in all of these halo poten- tials is perhaps unrealistically high, but chaotic orbits canbe found in more realistic potentials as well, as in Hunter2005.) It is worth noting that in this case, the thickeningof the stream arises within less than five radial oscillations.Long integration times are apparently not required to seethe influence of chaos on tidal streams.The results in Figure 11 demonstrate we can obtain rea-sonable results from the particle-spray model, even when theprogenitor orbit is irregular, though some areas of disagree-ment can be found in both rows. In fact, having first noticeda highly thickened stream in a similar orbit and potential us-ing the particle-spray model, we suspected a bug in our codeuntil the N -body model demonstrated the same behavior.For the spray model to work on a chaotic orbit, the parti-cle coordinates must be evolved by direct orbit integrationrather than by using an action-angle formalism.The fraction of orbits exhibiting chaotic behavior is gen-erally small, though nonzero, in the idealized axisymmetricmodels often used in stream modeling (Hunter 2005). How-ever, the fraction of orbits demonstrating chaotic behaviorin realistic galaxy potentials, which include triaxiality, ra-dial shape dependence, time dependence, and substructure,may well be higher. Also, real galaxy potentials and hencethe orbits within them evolve with time. In principle, or-bits may sweep through chaotic regions, puffing up their c (cid:13) , 1–20 M. A. Fardal et al. streams. This would raise the fraction of tidal streams af-fected by chaotic behavior beyond the fraction of orbits thatare currently chaotic.In this last set of simulations, we chose the parameterssuch that the satellite remains partly intact until the end ofthe simulation. If the satellite is more diffuse, it comes apartearlier in the N -body run than in the particle-spray model.Presumably this is due to disk shocking, though we havenot investigated the cause in detail. The particles are alsospread more evenly along the stream than in the particle-spray model, where there tends to be a central hole. Thesediscrepancies join those mentioned earlier as the main dis-crepancies with the N -body results: differences in the tailsof the particle distributions, inaccuracy in the satellite massas a function of time, and inaccuracy of the streams nearthe satellite Lagrange points. However, overall the agree-ment between these two methods is impressive, suggestingthe stream recipe is suitable for fast modeling of streams ina Bayesian context. We have presented a recipe for modeling tidal streams with-out the cost of a full N -body simulations. The typical com-putation times are a few minutes, as opposed to a typicalcomputation time of six hours for the N -body runs in Ta-ble 1. This recipe has many similarities with recent sugges-tions by Gibbons et al. (2014) and Bonaca et al. (2014), butdiffers in several respects. Here we discuss some aspects ofour recipe in the context of these and other methods.Our method clearly reproduces the differences betweenthe stream and the orbit that were reviewed in Section 2.2.As discussed in that section, the stream behavior is deter-mined primarily by the distribution of orbital frequencies(or equivalently the actions) upon release. Even though wehave made simplifying assumptions about the way the starsare released, and neglected the subsequent influence of thesatellite, the behavior of streams is still remarkably faithfulto the N -body results.Our model also reproduces the dispersion along vari-ous dimensions of the stream fairly well. Streakline tech-niques such as those of Varghese et al. (2011) or K¨upperet al. (2012) by design do not model this dispersion. There-fore, they fail in particular to reproduce the portions of thestream furthest from the satellite, which are dictated by theextremes of the frequency distribution at the first disruptivepericenter. They similarly fail to reproduce the extent of anyother lobes or streaky features within the stream. Thus whilethey are certainly easier to compute, the folded streaklinesevidenced by K¨upper et al. (2012) are not necessarily easierto fit to the stream than a full particle distribution. In ad-dition, the length and width of the stream contains usefulinformation on the satellite mass and overall duration of thedisruption process. These parameters are in turn valuable inobtaining a more accurate fit to the orbit and the potential.Substructure in streams provides important constraints,as described in section 2.3 and further illustrated in sec-tion 5. A crucial point in modeling these streams is that theactions and frequencies of released stars must vary over theradial cycle. In the recipe of Bovy (2014) for example, thedistribution of actions and frequencies is assumed constant. This type of model is useful in modeling old streams wherethe substructure overlaps, but is unable to reproduce thecomplex spatial pattern seen in for example Figure 10. Theejection recipe in our model naturally reproduces the varia-tion in orbital properties. In principle, this variation can beincorporated into action-angle stream models as well, whichwould be useful in obtaining large particle numbers (or accu-rate stream distribution functions) at a fixed computationalcost.While an action-angle approach has important advan-tages, irregular orbits will pose a problem for this approach.For these orbits, directly integration of the orbits as we havedone here seems the only viable method. The results shownin Figure 11 suggest that chaotic orbits result in puffed-upstreams within a few orbital times. Further investigation intothe effects of chaos on realistic galactic potentials and tidalstreams within them seems warranted.Our model has a varying emission rate peaking nearpericenter, in contrast to the constant emission used byBonaca et al. (2014) and Bovy (2014) among others. We havetested the effect of a uniform emission rate versus the sim-ulation in Figure 10. The effects are more subtle than thoseof the oscillating actions and frequencies; the substructureis still easily visible. However, the influence of the emissionrate is still evident in the abundance of particles betweenthe clearly separated streams, which fill in the distinct shellsseen at apocenter of the radial loops and blur the separatestreams visible in the energy plot. Thus this effect too isimportant in modeling streams with visible substructure.Our model also explicitly takes into account the de-crease of the satellite mass with time. This improves uponthe treatment in the models of Gibbons et al. (2014) andBonaca et al. (2014), as well as the earlier streakline mod-els, where the mass is assumed fixed. Models with multipleinteractions can be modeled fairly accurately with our sim-ple assumptions, as long as the satellites and orbits resemblethose used in our tests. For less disruptive encounters, morecareful modeling of the tidal excitation is probably neces-sary. Very disruptive encounters will require a criterion fortotal disruption and a prescription for the particles ejectedat that point. Including a disk in the host potential couldlead to disk shocking and a more rapid mass loss than spec-ified by our current recipe. Finally, different satellite pro-files will probably require a recalibration of the mass lossprescription. As discussed above, previous work has shownthat satellites with various density profiles evolve in a pre-dictable manner during moderate rates of tidal stripping, sothis recalibration should be possible as long as the satellitesare dynamically hot systems. The mass evolution in turnaffects the actions and frequencies of released particles, andignoring it in the case of Sagittarius for example may leadto significant biases.When compared to the Bonaca et al. (2014) method,two more differences stand out. First, they release stars atthe angular velocity of the satellite, which results in a largeroscillation of actions and frequencies over the radial cyclethan in our equations 6 and 20. Second, they blur the re-lease velocity of the satellite by the velocity dispersion in thesatellite, which in many cases is a larger dispersion than weuse. In some cases their assumptions work well, but in caseswith small or large eccentricities we find that they can pro-duce streams that are too smeared out in azimuth or radius. c (cid:13)000
63, the orbit be-comes regular again, and in tandem the stream becomesnarrow again so that it much more closely resembles the toppanel. (Admittedly, the flattening in all of these halo poten- tials is perhaps unrealistically high, but chaotic orbits canbe found in more realistic potentials as well, as in Hunter2005.) It is worth noting that in this case, the thickeningof the stream arises within less than five radial oscillations.Long integration times are apparently not required to seethe influence of chaos on tidal streams.The results in Figure 11 demonstrate we can obtain rea-sonable results from the particle-spray model, even when theprogenitor orbit is irregular, though some areas of disagree-ment can be found in both rows. In fact, having first noticeda highly thickened stream in a similar orbit and potential us-ing the particle-spray model, we suspected a bug in our codeuntil the N -body model demonstrated the same behavior.For the spray model to work on a chaotic orbit, the parti-cle coordinates must be evolved by direct orbit integrationrather than by using an action-angle formalism.The fraction of orbits exhibiting chaotic behavior is gen-erally small, though nonzero, in the idealized axisymmetricmodels often used in stream modeling (Hunter 2005). How-ever, the fraction of orbits demonstrating chaotic behaviorin realistic galaxy potentials, which include triaxiality, ra-dial shape dependence, time dependence, and substructure,may well be higher. Also, real galaxy potentials and hencethe orbits within them evolve with time. In principle, or-bits may sweep through chaotic regions, puffing up their c (cid:13) , 1–20 M. A. Fardal et al. streams. This would raise the fraction of tidal streams af-fected by chaotic behavior beyond the fraction of orbits thatare currently chaotic.In this last set of simulations, we chose the parameterssuch that the satellite remains partly intact until the end ofthe simulation. If the satellite is more diffuse, it comes apartearlier in the N -body run than in the particle-spray model.Presumably this is due to disk shocking, though we havenot investigated the cause in detail. The particles are alsospread more evenly along the stream than in the particle-spray model, where there tends to be a central hole. Thesediscrepancies join those mentioned earlier as the main dis-crepancies with the N -body results: differences in the tailsof the particle distributions, inaccuracy in the satellite massas a function of time, and inaccuracy of the streams nearthe satellite Lagrange points. However, overall the agree-ment between these two methods is impressive, suggestingthe stream recipe is suitable for fast modeling of streams ina Bayesian context. We have presented a recipe for modeling tidal streams with-out the cost of a full N -body simulations. The typical com-putation times are a few minutes, as opposed to a typicalcomputation time of six hours for the N -body runs in Ta-ble 1. This recipe has many similarities with recent sugges-tions by Gibbons et al. (2014) and Bonaca et al. (2014), butdiffers in several respects. Here we discuss some aspects ofour recipe in the context of these and other methods.Our method clearly reproduces the differences betweenthe stream and the orbit that were reviewed in Section 2.2.As discussed in that section, the stream behavior is deter-mined primarily by the distribution of orbital frequencies(or equivalently the actions) upon release. Even though wehave made simplifying assumptions about the way the starsare released, and neglected the subsequent influence of thesatellite, the behavior of streams is still remarkably faithfulto the N -body results.Our model also reproduces the dispersion along vari-ous dimensions of the stream fairly well. Streakline tech-niques such as those of Varghese et al. (2011) or K¨upperet al. (2012) by design do not model this dispersion. There-fore, they fail in particular to reproduce the portions of thestream furthest from the satellite, which are dictated by theextremes of the frequency distribution at the first disruptivepericenter. They similarly fail to reproduce the extent of anyother lobes or streaky features within the stream. Thus whilethey are certainly easier to compute, the folded streaklinesevidenced by K¨upper et al. (2012) are not necessarily easierto fit to the stream than a full particle distribution. In ad-dition, the length and width of the stream contains usefulinformation on the satellite mass and overall duration of thedisruption process. These parameters are in turn valuable inobtaining a more accurate fit to the orbit and the potential.Substructure in streams provides important constraints,as described in section 2.3 and further illustrated in sec-tion 5. A crucial point in modeling these streams is that theactions and frequencies of released stars must vary over theradial cycle. In the recipe of Bovy (2014) for example, thedistribution of actions and frequencies is assumed constant. This type of model is useful in modeling old streams wherethe substructure overlaps, but is unable to reproduce thecomplex spatial pattern seen in for example Figure 10. Theejection recipe in our model naturally reproduces the varia-tion in orbital properties. In principle, this variation can beincorporated into action-angle stream models as well, whichwould be useful in obtaining large particle numbers (or accu-rate stream distribution functions) at a fixed computationalcost.While an action-angle approach has important advan-tages, irregular orbits will pose a problem for this approach.For these orbits, directly integration of the orbits as we havedone here seems the only viable method. The results shownin Figure 11 suggest that chaotic orbits result in puffed-upstreams within a few orbital times. Further investigation intothe effects of chaos on realistic galactic potentials and tidalstreams within them seems warranted.Our model has a varying emission rate peaking nearpericenter, in contrast to the constant emission used byBonaca et al. (2014) and Bovy (2014) among others. We havetested the effect of a uniform emission rate versus the sim-ulation in Figure 10. The effects are more subtle than thoseof the oscillating actions and frequencies; the substructureis still easily visible. However, the influence of the emissionrate is still evident in the abundance of particles betweenthe clearly separated streams, which fill in the distinct shellsseen at apocenter of the radial loops and blur the separatestreams visible in the energy plot. Thus this effect too isimportant in modeling streams with visible substructure.Our model also explicitly takes into account the de-crease of the satellite mass with time. This improves uponthe treatment in the models of Gibbons et al. (2014) andBonaca et al. (2014), as well as the earlier streakline mod-els, where the mass is assumed fixed. Models with multipleinteractions can be modeled fairly accurately with our sim-ple assumptions, as long as the satellites and orbits resemblethose used in our tests. For less disruptive encounters, morecareful modeling of the tidal excitation is probably neces-sary. Very disruptive encounters will require a criterion fortotal disruption and a prescription for the particles ejectedat that point. Including a disk in the host potential couldlead to disk shocking and a more rapid mass loss than spec-ified by our current recipe. Finally, different satellite pro-files will probably require a recalibration of the mass lossprescription. As discussed above, previous work has shownthat satellites with various density profiles evolve in a pre-dictable manner during moderate rates of tidal stripping, sothis recalibration should be possible as long as the satellitesare dynamically hot systems. The mass evolution in turnaffects the actions and frequencies of released particles, andignoring it in the case of Sagittarius for example may leadto significant biases.When compared to the Bonaca et al. (2014) method,two more differences stand out. First, they release stars atthe angular velocity of the satellite, which results in a largeroscillation of actions and frequencies over the radial cyclethan in our equations 6 and 20. Second, they blur the re-lease velocity of the satellite by the velocity dispersion in thesatellite, which in many cases is a larger dispersion than weuse. In some cases their assumptions work well, but in caseswith small or large eccentricities we find that they can pro-duce streams that are too smeared out in azimuth or radius. c (cid:13)000 , 1–20 eneration of mock tidal streams The streams in Bonaca et al. (2014) are both produced andanalyzed with the same method, so it is unclear whetherthese assumptions have much effect on the larger conclu-sions in their paper. Still, these differences are only stronglyevident in certain situations. Overall, our results stronglyvalidate their general approach of modeling particle spraysusing orbits in the host potential alone.The method of Gibbons et al. (2014) is somewhat dif-ferent, in that they include an ongoing force from the satel-lite in calculating the particle evolution. It thus closely re-sembles a standard restricted N -body method. (A simplealternative to their method is to run a restricted N -bodysimulation after first removing the central particles of theprogenitor, because they remain bound to the satellite andadd nothing but computation time.) In their runs, betweena quarter and a half of the released particles are recapturedby the satellite and later released in bursts at pericentricpassage. Our method has the advantage that we do notneed to calculate the evolution of these recaptured parti-cles, which require short timesteps compared to the freelymoving stream particles. Gibbons et al. (2014) find that theforce from the satellite must be included. However, we findthat it can be omitted with an appropriate choice of releaseconstants (equations 19–24); the time-integrated effect of theprogenitor is encapsulated by our larger action offset. Onecan then use a simpler, general-use package for calculatingthe orbital evolution, rather than more specialized codes. Inour case we have used the galpy package which contains alarge assortment of potentials when orbit integrations arerequired.Another difference is that our method prescribes therate of particle release as a function of orbital phase. Whilethis form may not be completely accurate, it is unclearwhether the alternative of using recaptured particles to gen-erate bursts biases the results in some way. Finally, the ve-locity dispersion of released particles is prescribed by ourrecipe, and tied to the satellite mass and other parametersof the problem, where in Gibbons et al. (2014) it is taken tobe a free fitting parameter. Our method thus discards lessinformation, in principle increasing the statistical power ofthe fit. The main advantage of the Gibbons et al. (2014)method is the automatic inclusion of multiple interactionsbetween the satellite and its stream, which can create streamgaps like those visible in Figure 10. Also, at very large satel-lite masses the effect of the satellite will not converge withinan orbital wrap and this will alter the morphology of thetidal stream (Choi et al. 2007), making explicit inclusion ofthe satellite force necessary.Of course, one may question to what extent the ideal-ized N -body simulations used in this paper reflect the casesexpected in reality. Bonaca et al. (2014) describes an inter-esting experiment, where streams were generated in bothfixed potentials and in a cosmological simulation of a halo,and then fitted to constrain a parameterized model of thedark matter halo. While the fixed potential streams con-sistently returned accurate results, the cosmological halostreams yielded highly biased or imprecise results in manycases. Bonaca et al. suggest the differences arise from grad-ual evolution in the potential, interaction with subhalos,and deviations of the true potential shape from the fittedform. We also note the effect of chaotic orbits, as in Fig-ure 11 above, may play a role in these differences. There are some reasons to think the results of Bonaca et al. aretoo pessimistic about the scientific yield from studying tidalstreams. Neither of the two “streams” displayed from theirmuch larger sample particularly resembles prominent MilkyWay streams in the cosmological case, while one of themdoes not look much like a stream in the fixed-potential caseeither. We suggest that in the outer halo, progenitors thatare more massive than the globular-like systems they em-ploy (mass of 2 × M (cid:12) ) are necessary in order to gener-ate streams that yield robust measures of the potential. (Thefixed-potential case does yield good results, but this presum-ably depends on the assumption of perfect six-dimensionalobservational data.) This would align with earlier work usingmore massive progenitors (Siegal-Gaskins & Valluri 2008),which found only minor effects of substructures on tidalstreams. Still, the analysis of Bonaca et al. (2014) raisesinteresting challenges to the project of stream fitting, andfurther work on how to interpret and combine results fromindividual stream fits is clearly required.To the challenges facing stream models, we may add theissue of dynamical friction. Stream models are often initial-ized with a satellite at a particular place and time, usuallya particular pericentric passage. Dynamical friction is fre-quently invoked as a way of bringing this satellite onto itsdestructive orbit, but it is then often neglected in the actualcalculation of the stream, or else applied only to the satellitebut not the stream. Full N -body simulations with a live hostcan implement dynamical friction accurately. However, thistechnique is extremely expensive, since high particle num-bers are necessary to suppress the effects of particle noise.Thus at this point, a good technique for incorporating dy-namical friction consistently in stream models has yet to beintroduced. In summary, we have presented and tested a recipe for mod-eling tidal streams as a collection of particles, released andevolved in the host potential without the added influenceof the progenitor satellite. This recipe includes a prescrip-tion for calculating the mass lost from the progenitor, andthus the mass in the tidal stream as well. The mass lossand orbital properties of the debris both oscillate with thesatellite position, producing substructure within the tidalstream. The recipe assumes hot, compact one-componentprogenitors, and will require adaptation for more compli-cated cases, or orbits very different from the moderatelydisrupting cases examined here.The intended use of this recipe is for Bayesian sam-pling, which requires a very large number of trial modelsto obtain accurate constraints on parameters. The particle-spray method adds one more arrow to the quiver for thoseseeking to use tidal streams to obtain physical insights. Theobservational situation has been developing rapidly. Numer-ous streams have been found in the Milky Way using theSloan survey, and in M31 with the PAndAS survey. Newtechniques have revealed tidal streams in more distant galax-ies with even modest-sized telescopes, and kinematic tracershave been demonstrated in some of these streams as well.Forthcoming results from Pan-STARRS and Gaia, amongother surveys, should transform the field of tidal streams into c (cid:13) , 1–20 M. A. Fardal et al. a high-precision industry. We expect particle-spray modelingsuch as that tested here to be one of its chief tools.
ACKNOWLEDGMENTS
We thank Tom Quinn and Joachim Stadel for the use of
PKDGRAV , and Josh Barnes for the use of
ZENO . We thankJo Bovy for making public the galpy package and assistingwith its use, We also thank Jo Bovy and Aaron Romanowskyfor interesting discussions of tidal stream physics and ob-servations. MAF, MDW, and SH acknowledge support byNSF grant AST-1009652 to the University of Massachusetts.MDW also acknowledges the support of NSF AST-0907951.
REFERENCES
Aarseth S. J., Heggie D. C., 1998, MNRAS, 297, 794Belokurov V., Zucker D. B., Evans N. W., Gilmore G., OtherS., Authors D. M., 2006, ApJ, 642, L137Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edi-tion. Princeton University PressBonaca A., Geha M., Kuepper A. H. W., Diemand J., JohnstonK. V., Hogg D. W., 2014, ArXiv e-printsBovy J., 2014, ArXiv e-printsCarlberg R. G., 2013, ApJ, 775, 90Chernoff D. F., Weinberg M. D., 1990, ApJ, 351, 121Choi J.-H., Weinberg M. D., Katz N., 2007, MNRAS, 381, 987Eyre A., Binney J., 2011, MNRAS, 413, 1852Fardal M. A., Babul A., Geehan J. J., Guhathakurta P., 2006,MNRAS, 366, 1012Fardal M. A., Weinberg M. D., Babul A., Irwin M. J.,Guhathakurta P., Gilbert K. M., Ferguson A. M. N., IbataR. A., Lewis G. F., Tanvir N. R., Huxor A. P., 2013, MNRAS,434, 2779Gibbons S. L. J., Belokurov V., Evans N. W., 2014, ArXiv e-printsGieles M., Alexander P. E. R., Lamers H. J. G. L. M., Baum-gardt H., 2014, MNRAS, 437, 916Gnedin O. Y., Ostriker J. P., 1997, ApJ, 474, 223Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn T.,2003, ApJ, 584, 541Heggie D. C., 2001a, in Steves B. A., Maciejewski A. J., eds,The Restless Universe Escape in Hill’s problem. pp 109–128Heggie D. C., 2001b, in Deiters S., Fuchs B., Just A., SpurzemR., Wielen R., eds, Dynamics of Star Clusters and the MilkyWay Vol. 228 of Astronomical Society of the Pacific ConferenceSeries, Mass Loss from Globular Clusters. p. 29Helmi A., White S. D. M., 1999, MNRAS, 307, 495Howley K. M., Geha M., Guhathakurta P., Montgomery R. M.,Laughlin G., Johnston K. V., 2008, ApJ, 683, 722Hunter C., 2005, Annals of the New York Academy of Sciences,1045, 120Johnston K. V., 1998, ApJ, 495, 297Just A., Berczik P., Petrov M. I., Ernst A., 2009, MNRAS, 392,969King I. R., 1966, AJ, 71, 64K¨upper A. H. W., Kroupa P., Baumgardt H., Heggie D. C., 2010,MNRAS, 401, 105K¨upper A. H. W., Lane R. R., Heggie D. C., 2012, MNRAS,420, 2700K¨upper A. H. W., MacLeod A., Heggie D. C., 2008, MNRAS,387, 1248Lamers H. J. G. L. M., Baumgardt H., Gieles M., 2010, MNRAS,409, 305Murali C., Weinberg M. D., 1997, MNRAS, 291, 717 Odenkirchen M., Grebel E. K., Dehnen W., Rix H., YannyB., Newberg H. J., Rockosi C. M., Mart´ınez-Delgado D.,Brinkmann J., Pier J. R., 2003, AJ, 126, 2385Pe˜narrubia J., Navarro J. F., McConnachie A. W., 2008, ApJ,673, 226Price-Whelan A. M., Hogg D. W., Johnston K. V., Hendel D.,2014, ArXiv e-printsSanders J. L., 2014, MNRAS, 443, 423Sanders J. L., Binney J., 2013a, MNRAS, 433, 1813Sanders J. L., Binney J., 2013b, MNRAS, 433, 1826Sanderson R. E., Helmi A., Hogg D. W., 2014, in Feltzing S.,Zhao G., Walton N. A., Whitelock P., eds, IAU SymposiumVol. 298 of IAU Symposium, Action-space clustering of tidalstreams to map the Galactic potential. pp 207–212Siegal-Gaskins J. M., Valluri M., 2008, ApJ, 681, 40Sohn S. T., van der Marel R. P., Carlin J. L., Majewski S. R.,Kallivayalil N., Law D. R., Anderson J., Siegel M. H., 2014,ArXiv e-printsTakahashi K., Portegies Zwart S. F., 2000, ApJ, 535, 759van den Bosch F. C., Lewis G. F., Lake G., Stadel J., 1999, ApJ,515, 50Varghese A., Ibata R., Lewis G. F., 2011, MNRAS, 417, 198c (cid:13)000