Predicted Extension of the Sagittarius Stream to the Milky Way Virial Radius

Galaxy formation is a continuous process, and as a result, the formation redshift of galaxies such as the MW and Sgr dwarf cannot be pinpointed exactly. We aim to trace the system’s history out to a redshift of z = 1, corresponding to approximately 8 Gyr of lookback time. This choice is based on the age of the M-giants in the stream, estimated by Bellazzini et al. (2006) to be 8.0 ± 1.5 Gyr. We adopt this evolutionary timescale as the guiding principle for our initial conditions, postulating that the Sgr progenitor may first have crossed the MW’s virial radius around that time. We adopt a fiducial virial mass of 10¹² M_⊙ for the MW, consistent with recent estimates (e.g., Xue et al. 2008). Various simulations of halo assembly histories suggest that today’s MW-like haloes would have built up at least half of their mass by z = 1, with a significant scatter of ∼2 × 10¹¹ M_⊙ depending on the halo’s accretion history (e.g., Torrey et al. 2015; Lu et al. 2016). Using the spherically symmetric top-hat framework of halo formation, the virial radius of a 5 × 10¹¹ M_⊙ galaxy collapsing at z = 1 is approximately 124 kpc (Barkana & Loeb 2001). We therefore initiate the Sgr progenitor at a Galactocentric radius of 125 kpc and aim to trace its evolution over a period of 8 Gyr.

Halo growth is believed to occur inside-out, with later additions of mass being appended on the outskirts of the halo (e.g., Loeb & Peebles 2003). This picture has been validated in simulations (e.g., Wellons et al. 2016) and observations (e.g., de la Rosa et al. 2016) of high-redshift massive compact galaxies. Cosmological simulations also show that after z = 1, most MW analogs do not undergo a major merger (e.g., Fakhouri et al. 2010) and have inner halo profiles that remain essentially fixed (e.g., Gao et al. 2004; Mollitor et al. 2015). Studies of the stellar populations in the MW disk suggest a quiet evolution since z = 2, with no significant mergers in the last ∼10 Gyr (Wyse 2001; Hammer et al. 2007). The Sgr orbit in our model is contained inside Galactocentric radii <125 kpc and is consequently not impacted by the later addition of matter outside this sphere (assuming spherical symmetry). As a result, we maintain the same initial halo profile throughout each simulation. A halo’s scale radius r_s is related to its virial radius R_vir and concentration parameter c by the relationship r_s = R_vir/c. As the halo accretes mass on its outskirts, the scale radius is maintained constant while the virial radius and the concentration parameter grow. Following the growth in concentration of galaxies in the models of Diemer & Kravtsov (2015), an MW analog would have grown from c ≃ 7, R_vir ≃ 125 kpc at z = 1 to c ≃ 10, R_vir ≃ 200 kpc at z = 0, giving a nearly constant scale radius of ∼18–20 kpc. For a z = 0 MW analog (with parameters as outlined in Table 1), the mass inside a radius of 125 kpc is approximately 7 × 10¹¹ M_⊙, consistent with the virial mass of (5 ± 2) × 10¹¹ M_⊙ of its z = 1 progenitor. Beyond the prescriptions for the MW’s evolution adopted here, Peñarrubia et al. (2006) provide an in-depth study of modeling tidal streams in evolving dark matter haloes. Their simulations suggest that the debris configuration reflects the present galactic potential shape rather than its evolution. Our choice of a constant host potential profile throughout the orbit of Sgr is therefore unlikely to have a strong effect on the phase-space distribution of the stream.

The MW potential includes three components: a dark matter halo and bulge following Hernquist profiles (Hernquist 1990), and an exponential disk:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{MW}}={{\rm{\Phi }}}_{\mathrm{halo}}+{{\rm{\Phi }}}_{\mathrm{disk}}+{{\rm{\Phi }}}_{\mathrm{bulge}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{halo}}(r)=-\displaystyle \frac{{{GM}}_{\mathrm{halo}}}{r+{r}_{{\rm{H}}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{disk}}(r)=-{{GM}}_{\mathrm{disk}}\left(\displaystyle \frac{1-{e}^{-r/{b}_{0}}}{r}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{bulge}}(r)=-\displaystyle \frac{{{GM}}_{\mathrm{bulge}}}{r+{c}_{0}}.\end{eqnarray} \tag{ 4 }$$

For simplicity, the disk and bulge components are omitted for Sgr in the semianalytic model, but included in the full simulation described in Section 3. The parameters of both galactic potentials are summarized in Table 1. The parameters of the Hernquist halo potentials are chosen such that the mass enclosed within the radius of interest matches that of fiducial Navarro, Frenk, & White (NFW; Navarro et al. 1997) profiles. The resulting curves of enclosed mass are shown in Figure 1. For the MW, we adjust the enclosed mass inside the initial distance between the galaxy centers, d_init = 125 kpc. A total halo mass of M_halo = 1.25 × 10¹² M_⊙ and a scale radius of r_H = 38.35 kpc match the fiducial NFW halo with a virial mass M_MW = 10¹² M_⊙ and a concentration parameter of 10. For Sgr, we use M_halo = 1.3 × 10¹⁰ M_⊙ and r_H = 9.81 kpc. This corresponds to the mass enclosed within the initial tidal radius, ${r}_{\mathrm{tidal},\mathrm{init}}\simeq 25\,\mathrm{kpc}$, of an NFW potential with M_Sgr = 10¹⁰ M_⊙ and c_Sgr = 8.

Zoom In Zoom Out Reset image size

Our choice of initial parameters for Sgr is based on the study by Niederste-Ostholt et al. (2010), who reconstruct the properties of the progenitor by conducting a census of the stellar tidal debris. A lower bound of (9.6–13.2) × 10⁷ L_⊙ is inferred for the progenitor’s luminosity by summing up the luminosities of the Sgr core, and leading and trailing arms. Relating this value to results from cosmological N-body simulations, Niederste-Ostholt et al. (2010) estimate a mass of ∼10¹⁰ M_⊙ for the Sgr dark matter halo prior to tidal disruption. Based on this choice of progenitor mass, the models of the mass–concentration relation of Diemer & Kravtsov (2015) suggest a low concentration of c ≃ 8 at z = 1. We then use these parameters as our reference NFW profile, and tune the Hernquist profile parameters to match the mass enclosed inside the initial tidal radius of the system (see Figure 1). The resulting Hernquist scale radius r_H = 9.81 kpc is consistent with the values used by Gómez et al. (2015) for their “Light Sgr” model, which has a mass of 3.2 × 10¹⁰ M_⊙. Since only one wrap of the debris is considered and additional luminous matter may be located at apocentric pile-ups, the true Sgr progenitor mass may exceed the lower bound established by Niederste-Ostholt et al. (2010). Given a total luminosity of ∼10⁸ L_⊙, cosmological abundance matching would suggest a higher mass of ∼10^10.5–10¹¹ M_⊙ (Conroy & Wechsler 2009; Behroozi et al. 2010). Exploring a range of different progenitor models, Gibbons et al. (2017) find that masses ≳6 × 10¹⁰ M_⊙ are most consistent with the velocity dispersion of the stellar populations in the stream. Still, the orbit presented here is the first physically motivated model for the Sgr infall that includes a progenitor mass closer to the most recent mass estimates.

In our search for an orbital model for the Sgr dwarf spheroidal, we start testing possible trajectories with an approximate point-particle approach. From predetermined initial conditions, the equations of motion are solved numerically forward in time with a fourth-order Runge–Kutta method. Both the dynamics of Sgr in the MW potential and the MW’s response to the gravitational attraction of Sgr are modeled, including tidal stripping and dynamical friction for Sgr. We emphasize the fact that the MW potential is not held fixed at the origin; rather, we account for the mutual gravitational attraction between the MW and Sgr and allow the MW to move freely about the common center of mass of the system. As the studies of Dierickx et al. (2014) and Gómez et al. (2015) have shown, the response of an MW-like host galaxy to an infalling satellite can be significant. This is particularly relevant in the regime of high separation and high Sgr progenitor mass explored here.

Keeping the initial separation fixed, the remaining free parameters quantify the amount of initial orbital angular momentum we grant Sgr at its starting location. We specify two quantities: v_init, the magnitude of the Sgr velocity, and θ_init, the angle between the velocity vector and the direction to the MW center (see Figure 2). We integrate a two-dimensional grid of 2500 orbits with v_init ranging from 0 to ∼320 km s⁻¹, the NFW escape velocity from the MW halo at the initial radius of the Sgr progenitor, and θ_init ranging from 10° to 90°. With total energy per unit mass defined as $E=\tfrac{1}{2}{v}_{\mathrm{init}}^{2}-{\rm{\Phi }}({d}_{\mathrm{init}})$, the upper limit on velocity investigated here corresponds to E = 0, i.e., a marginally bound Sgr falling into the MW from a larger distance. We specify a nonzero lower limit on the initial angle in order to avoid purely radial orbits, which cause pathological behavior in the model. Angles >90° are excluded as they would imply an Sgr progenitor velocity directed away from the MW.

Zoom In Zoom Out Reset image size

Throughout the orbital integration, the MW experiences acceleration from a standard Hernquist gravitational potential for a tidally truncated Sgr. For the Sgr satellite, however, the following acceleration is calculated at Galactocentric coordinates ${\boldsymbol{r}}$:

$$\begin{eqnarray}&&{\boldsymbol{a}}=-\displaystyle \frac{{{GM}}_{\mathrm{halo}}(\lt r)}{{r}^{3}}{\boldsymbol{r}}-{\rm{\nabla }}{{\rm{\Phi }}}_{\mathrm{disk}}-{\rm{\nabla }}{{\rm{\Phi }}}_{\mathrm{bulge}}+{{\boldsymbol{a}}}_{\mathrm{DF}},\end{eqnarray} \tag{ 5 }$$

where M_halo(<r) is the mass interior to radius r of the host halo, given simply by ${M}_{\mathrm{halo}}{(\lt r)={M}_{\mathrm{halo}}[{r}^{2}/(r+{r}_{H})}^{2}]$ (Hernquist 1990). The acceleration introduced by dynamical friction, ${{\boldsymbol{a}}}_{\mathrm{DF}}$, is modeled using the Chandrasekhar formula (Binney & Tremaine 1987, p. 747, Equation (7.18)):

$$\begin{eqnarray}{{\boldsymbol{a}}}_{\mathrm{DF}} & = & -\displaystyle \frac{4\pi \mathrm{ln}({\rm{\Lambda }}){G}^{2}\rho (r){M}_{\mathrm{Sgr}}(\lt {r}_{t})}{{v}^{3}}\\ & & \times \left[\mathrm{erf}(X)-\displaystyle \frac{2X}{\sqrt{\pi }}\exp (-{X}^{2})\right]{\boldsymbol{v}},\end{eqnarray} \tag{ 6 }$$

where ρ(r) is the density of the host halo at Galactocentric distance r, M_Sgr(<r_t) is the Sgr mass interior to its minimum tidal radius r_t, and $\mathrm{ln}({\rm{\Lambda }})$ is the Coulomb logarithm. The remaining term involves $X=v/\sqrt{2}\sigma$, where v is the satellite velocity and σ is the one-dimensional velocity dispersion of particles in the host halo (given by Hernquist 1990, Equation (10)). Based on the formalism used by Besla et al. (2007) in their study of the orbital evolution of the Magellanic Clouds, we adopt an alternative time-dependent Coulomb logarithm as follows (Hashimoto et al. 2003):

$$\begin{eqnarray}&&\mathrm{ln}({\rm{\Lambda }})=\mathrm{ln}\left(\displaystyle \frac{r}{1.4\epsilon }\right),\end{eqnarray} \tag{ 7 }$$

where is a softening length variable used by Hashimoto et al. (2003) to model the Large Magellanic Cloud (LMC) with a Plummer sphere. This time-dependent parameterization of the Coulomb logarithm is also in agreement with the calibration carried out by van der Marel et al. (2012a) for the M31–M33 interaction (which is of similarly unequal mass). We ran a series of semianalytic orbits to test the friction parameterization, and find that setting  ≃ 1 kpc gives the best agreement with N-body integrations of the same initial conditions. This halo parameter is smaller by a factor of 3 than the parameterization used by Hashimoto et al. (2003) for the LMC, which is not surprising given that the mass ratio between host and satellite is more unequal by a factor of 10 in our case. Comparing to the standard formalism where Λ = b_max/b_min (Binney & Tremaine 1987, p. 747), we have in essence set the impact parameter b_min to ∼1.4 kpc (a reasonable value given the large range of initial conditions explored here) and introduced a time dependence for the cutoff radius b_max. These adjustments are generally expected to have a minor effect given their logarithmic contribution to the dynamical friction exerted on Sgr.

Sgr experiences tidal forces as it moves across the MW halo, leading to the outermost layers of material being stripped. The instantaneous tidal radius for Sgr, r_t, is found by solving the following equation numerically (King 1962):

$$\begin{eqnarray}&&{r}_{t}=r{\left[\displaystyle \frac{1}{2}\displaystyle \frac{{M}_{\mathrm{Sgr}}(\lt {r}_{t})}{{M}_{\mathrm{MW}}(\lt r)}\right]}^{1/3}.\end{eqnarray} \tag{ 8 }$$

Since material that has been stripped is considered lost, the tidal radius is computed at each time step and, if necessary, updated to the lowest value calculated so far. Using this prescription, we model dynamical friction on a progressively less massive Sgr halo. Since the lost Sgr mass is negligible compared to the MW mass, we ignore its contribution to the MW potential.

From its initial location at d_init = 125 kpc with varying velocity vector ${{\boldsymbol{v}}}_{\mathrm{init}}$, the orbit of the Sgr progenitor is integrated forward in time for 10 Gyr with 10 Myr time steps. The current position of Sgr can be described by its Galactic longitude and latitude (l, b) = (56, −142) (Majewski et al. 2003) and its heliocentric distance d_helio = 25 ± 2 kpc (Kunder & Chaboyer 2009). We define a right-handed Cartesian coordinate system centered on the Galactic Center, with the z-axis pointing toward the north Galactic pole and with the Sun’s current location at [−8, 0, 0] kpc (Honma et al. 2012; Reid et al. 2014). The coordinates of Sgr in this system can then be found with the following simple transformations (in units of kiloparsecs):

$$\begin{eqnarray}&&x={d}_{\mathrm{helio}}\cos b\cos l-8,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&y={d}_{\mathrm{helio}}\cos b\sin l,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&z=d\sin b.\end{eqnarray} \tag{ 11 }$$

If we apply these equations then the current position vector of the center of Sgr is ${{\boldsymbol{r}}}_{\mathrm{Sgr},\mathrm{obs}}=(16.1,2.35,-6.12)$ kpc, in agreement with the values provided, e.g., by Law & Majewski (2010). The heliocentric radial velocity of Sgr has been measured as 140 ± 0.33 km s⁻¹ (weighted mean of the estimates of the average velocities of Sgr,N and M54 in Table 5 of Bellazzini et al. 2008), and its proper motion in the equatorial coordinate system is (μ_α, μ_δ) = (−2.95 ± 0.18, −1.19 ± 0.16) mas yr⁻¹ (Massari et al. 2013). Using these heliocentric velocity components, the Cartesian, right-handed Galactic space velocity (U, V, W) is calculated following Johnson & Soderblom (1987). Finally, these heliocentric space velocities are converted to the Galactic Rest Frame (GSR) by adding contributions from the local standard of rest and solar peculiar motions (Schönrich et al. 2010; Reid et al. 2014):

$$\begin{eqnarray}{{\boldsymbol{v}}}_{\mathrm{GSR}}[\mathrm{km}\,{{\rm{s}}}^{-1}] & = & (U,V,W)+(0,237,0)\\ & & +(11.1,12.24,7.25),\end{eqnarray} \tag{ 12 }$$

This yields a current GSR velocity vector for Sgr of ${{\boldsymbol{v}}}_{\mathrm{Sgr},\mathrm{obs}}=(242.5,5.6,228.1)$ km s⁻¹ and a total magnitude of the velocity of 333 ± 30 km s⁻¹.

At every time step along the trajectory, the quality of the match to the present-day configuration is investigated. This is first done in a spherically symmetric sense before finding a best-match orientation for the Sun’s location. At every step, the MW–Sgr distance and the magnitude of the Sgr velocity in the GSR are compared to the observed values outlined above. For the times when Sgr is within 3σ of the correct Galactocentric distance and the correct magnitude of the velocity, the analysis progresses to the next step. At this stage, it is necessary to orient the system in order to further quantify the match to observables. We utilize the spherical symmetry of the host potential to identify which point on a sphere of 8 kpc radius around the Galactic Center provides the best match to the Sun’s location.

To this end, we construct a Fibonacci lattice of 1600 points with radius 8 kpc centered at the Galactic Center. The Fibonacci lattice is a convenient method to generate any number N of evenly distributed points on the surface of a sphere, with each point representing almost the same area (Gonzàlez 2010). Points are arranged along a tightly wound generative spiral, and use of the golden angle (∼1375) as the value for the longitudinal turn between consecutive points maximizes packing. With N = 1600, every lattice point occupies on average ∼25 deg², representing a one-dimensional positional uncertainty of approximately 5°. Scaled by the ratio of the Sgr heliocentric distance to the Sun’s Galactocentric distance (approximately a factor of 3), a lattice with N = 1600 points yields an angular uncertainty of ∼15–2° in the position of Sgr in the sky without becoming computationally prohibitive. This is acceptable given the many approximations and the uncertainties associated with the parameters of the galactic potentials.

For every candidate position of the Sun on the lattice, we compute the following six quantities: (i) the Galactic Center–Sgr distance; (ii) Sgr heliocentric distance; (iii) the Galactic Center–Sgr–Sun angle; (iv) the Sgr heliocentric radial velocity; (v) the magnitude of the Sgr heliocentric transverse velocity; (vi) the angle between the Sgr transverse velocity vector and the direction to the Galactic Center. The match to the corresponding measured numbers is quantified via a chi-squared test. The parameter values and associated errors used in the analysis are provided in Table 2. The lattice point with the lowest reduced chi-squared is identified and recorded as the best-match position for the Sun at each favorable time step in the run. The moment along each trajectory with the lowest reduced chi-squared between 7 and 8.5 Gyr after initialization is extracted in order to compare different sets of initial conditions.

Table 2.  Observable and Simulated Parameters Used in Chi-squared Analysis

Quantity	Observed Value	Semianalytic Model	N-body Model	Observational Error (σ)
(i) [kpc]	17.4	18.6	21.0	2.0
(ii) [kpc]	25.0	26.0	27.6	2.0
(iii) [deg]	6.92	7.6	10.9	5
(iv) [km s−1]	178.8	180.8	172.8	1.5
(v) [km s−1]	281.0	353.8	285.0	30
(vi) [deg]	162.0	160.4	162.8	1.5
Resulting ${\chi }_{r}^{2}$	⋯	1.32	2.28	⋯

Note. The six reference quantities used in the χ² analysis are provided along with the simulated parameters for both the best-match analytic and N-body models (see Section 3). The table rows are as follows: (i) Galactic Center (GC)–Sgr distance [kpc]; (ii) Sgr–Sun distance [kpc]; (iii) GC–Sgr–Sun angle [deg]; (iv) Sgr heliocentric line-of-sight velocity [km s⁻¹]; (v) Sgr heliocentric transverse velocity [km s⁻¹]; (vi) angle between Sgr transverse velocity vector and the direction to the GC [deg].

Download table as:  ASCIITypeset image

Figure 3 shows the reduced chi-squared values associated with the best-match snapshots along each trajectory, over the initial parameter space outlined in Section 2.2. In this color map, darker pixels indicate a closer match to the observed position and velocity of Sgr today. No matches are found for initial angles θ_init ≲ 30°, indicating that nearly radial orbits are incompatible with the distance and velocity constraints imposed by the available data for Sgr. At high initial incidence angles (θ_init ≳ 60°), a specific range of initial velocities (80 km s⁻¹ ≲ v_init ≲ 100 km s⁻¹) is preferred independently of the angle. Tangential velocities in this range are reasonable given analogous measurements of Local Group satellites: the M33 tangential velocity with respect to M31 is ∼129 km s⁻¹ (van der Marel et al. 2012b); the tangential velocities of the LMC and SMC with respect to the MW are ∼314 and ∼61 km s⁻¹, respectively (Kallivayalil et al. 2013).

Zoom In Zoom Out Reset image size

Contours of constant initial angular momentum are superimposed on the color map in Figure 3. Values of reduced chi-squared indicating a good match $({\mathrm{log}}_{10}({\chi }_{r}^{2})\lesssim 0.5)$ are tightly grouped together in a “valley” of constant angular momentum. This suggests that in order to reach the current observed position and velocity of Sgr, a narrow range of initial angular momenta is preferred. The v_init width of this range per unit mass and distance is approximately 20 km s⁻¹. We thus conclude that the model favors a specific value for the angular momentum Sgr had 7–8 Gyr ago, upon crossing the MW virial radius for the first time.

Next we extract the set of initial parameters with the lowest reduced chi-squared from Figure 3, and the resulting orbit is shown in Figure 4. In terms of the parameters defined in Figure 2, an initial velocity of v_init = 72.6 km s⁻¹ directed at an angle of θ_init = 808 yields a ${\chi }_{r}^{2}$ value of 1.32 at a time 7.71 Gyr after the start of the integration (see Table 2). The left panel of Figure 4 shows the resulting trajectory of the Sgr progenitor (in blue) and the MW’s drift (in black). The right panel shows the separation between the two galaxies as a function of time (dashed line) as well as the minimum tidal radius of the Sgr galaxy calculated as described in Section 2 (black line). Figure 5 presents a comparison of the best-match simulated and measured phase-space coordinates of Sgr. Overall the semianalytic model reproduces the position and velocity vector of Sgr very well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With the simultaneous aims of checking the semianalytic formalism described earlier and generating a mock Sgr stream, we utilize the N-body code GADGET (Springel 2005) to re-run the best-fit trajectory presented in Section 2.3. After initially exploiting a simplified spherically symmetric model for the host potential, we now take advantage of having pinned down the best-fit location for the Sun to introduce a flattened disk potential for the MW. With Sgr initialized at (125, 0, 0) kpc and with ${{\boldsymbol{v}}}_{\mathrm{init}}\simeq (-10,0,70)$ km s⁻¹, the MW disk inclination is defined by the Sun’s position at ∼(7, −1, 3) kpc. The parameters of the Hernquist halo, bulge, and exponential disk components of the simulated host and satellite galaxies are summarized in Table 1. Since the purpose of this study is to shed light on the dynamical history of Sgr (rather than model its star formation history, for example), the simulation does not include hydrodynamics. The absence of disk gas components should not affect the dynamics of the system, because gas represents only a small percentage of the galaxies’ mass budgets.

Full N-body simulations of the MW–Sgr interaction are costly because of the uneven ratio in progenitor mass between the two galaxies. As a result, only a few studies so far (e.g., Purcell et al. 2011, 2012; Gómez et al. 2015) have modeled the MW with a live halo rather than a static potential. However, given the significant mass of the Sgr satellite and its repeated passages at low Galactocentric distances, it is expected not only to have strong effects on the structure of the MW disk, but also to lead to significant overdensities of dark matter in the halo (Purcell et al. 2012). Both the host and satellite haloes are live in our simulation in order to capture such time-dependent effects on the structure of the MW potential. We aim to resolve the total visible mass of Sgr (∼10⁹ M_⊙—see Table 1) to the order of a few tens of thousands of particles in order to sufficiently populate the tidal stream. Therefore, we choose a mass resolution of 4 × 10⁴ M_⊙ per stellar particle, yielding 20,800 stellar particles in Sgr. This in turn implies a required ∼2.3 × 10⁶ particles to simulate the MW’s stellar mass of ∼9.4 × 10¹⁰ M_⊙. We use a particle resolution for dark matter of 10⁶ M_⊙, giving ∼1.2 × 10⁶ and ∼1.2 × 10⁴ dark matter particles for the MW and Sgr galaxies respectively. The majority of the computational cost in the simulation therefore comes from adequately resolving the baryonic component of the Sgr dwarf. The option of modeling the MW with fewer particles of higher mass is undesirable, because it leads to rapid two-particle relaxation and introduces artificial disruption of the Sgr satellite (Jiang & Binney 2000). We use an adaptive time step of maximum length 10 Myr, and the softening lengths for the baryonic and dark matter particles are 41 pc and 214 pc, respectively.

Figure 6 presents a general overview of the orbit of Sgr computed here. The apocentric and pericentric distances decrease at each passage due to dynamical friction. We note that friction is more efficient in the live simulation, gradually shrinking the apocentric distances reached by the satellite as compared to the semianalytic calculation. Similarly, the discrepancy between the orbital periods in each model grows in time, such that after 7 Gyr of evolution, the timing of the fourth pericenter passage differs by ∼0.6 Gyr. In the N-body case, the time corresponding to the present-day configuration is reached after 7.13 Gyr. At that time, the match between the observed and simulated phase-space coordinates of Sgr is quantified by a reduced chi-squared value of 2.28, slightly worse than for the semianalytic model (see Table 2). Figure 7 shows both the observed and simulated three-dimensional position and velocity of the Sgr core at the present time in relation to the Sun. Across the many test simulations carried out for this work, it is generally the case that the simulated Sgr core is slightly more distant than observed. This challenge in reaching small Galactocentric radii likely arises from our choice of initial conditions, with the Sgr progenitor starting at larger distances than previously considered by, for example, Law & Majewski (2010) and Purcell et al. (2011). Despite the much larger initial separation chosen here, the MW–Sgr distance and relative velocity vectors are still matched within approximately 2σ. Stronger dynamical friction may help in reaching smaller Galactocentric separations, suggesting that the Sgr progenitor may have initially been more massive than considered in this study (M_Sgr = 10¹⁰ M_⊙). Testing the dependence of the parameter match on the model MW halo could also provide useful constraints on its properties, a possibility we plan to explore in more depth in a follow-up paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 presents the behavior of the Sgr progenitor’s collisionless components at the time of best match in the simulation. Following Belokurov et al. (2014), the coordinates of the stream particles are projected on the debris plane defined by the pole located at (l_GC, b_GC) = (275°, −14°). Prominent tidal features are clearly visible for both the dark matter and stellar particles. Shells and arcs corresponding to apocentric pile-ups appear in both the stellar and dark matter components. Interestingly, both the dark and visible particles form large-scale streams at Galactocentric radii well beyond the apocentric distances reached by the Sgr core. The fact that some of the Sgr stars end up outside the initial orbit can be understood by analogy with the energy redistribution that occurs in tidal disruption events (TDEs). In TDEs, half of the stellar mass gains enough energy to reach the escape speed, while the rest becomes more tightly bound (see, e.g., Rees 1988; Guillochon et al. 2016). The energy gained by the leading arm from the tidal force could lead to large apocentric distances. Our simulation features such extended tidal arms because the initial orbital radius of Sgr is more remote than in previous studies. We have labeled the three most prominent distant extensions N_W, N_E, and S, standing for the northwestern, northeastern, and southern branches, respectively. We expect the sharpness of these stellar structures to be affected by the initial distribution of stars in the Sgr model. For a disk-less progenitor, the streams are likely to be less well defined. Currently, the most distant tidal debris identified are part of the trailing tail of the Sgr stream, with apocentric distances of the order of 102.5 ± 2.5 kpc according to Belokurov et al. (2014). Figure 6 shows that our model predicts significant tidal features at distances >100 kpc.

Zoom In Zoom Out Reset image size

In Figure 9 we present a comparison between available data and the modeled stream features projected on the Sgr orbital plane as defined by Belokurov et al. (2014). In the right panel, black crosshairs show data from Sloan Digital Sky Survey (SDSS) Data Releases 5 and 8 as aggregated in Figure 10 of Belokurov et al. (2014). The size of the markers is indicative of the largest distance errors present in the data set. Our choice of initial conditions for the Sgr baryonic components likely affects the thickness of the simulated stream. We did not try to vary the velocity dispersion of the stellar particles to optimize for stream width.

Zoom In Zoom Out Reset image size

The SDSS detections of both the leading and trailing arms have obvious counterparts in the simulation. The leading branch of the simulated stream presents a small offset both from the observed data and from the model by Law & Majewski (2010) shown in Belokurov et al. (2014, Figure 10): the apocentric distance appears ∼10 kpc larger than the measured value of 47.8 ± 0.5 kpc (Belokurov et al. 2014), representing a ∼20% mismatch. The position angle of the apocenter location is also slightly different from the data. Belokurov et al. (2014) have argued that the precession angle of only 932 ± 35 measured from the SDSS data is an indication that the density of MW dark matter falls more quickly with radius than a logarithmic potential (with typical precession angles of 120°). While orbital energy and angular momentum also play a minor role, the fact that the simulated precession angle appears ∼10° larger than the measured value might suggest that a steeper model potential is needed for the MW.

However, the model’s most promising feature is that it successfully reproduces the distant trailing arm of the Sgr stream. Initially detected by Newberg et al. (2003), remote Sgr stellar debris in the northern hemisphere were later confirmed by Ruhland et al. (2011), Drake et al. (2013), Belokurov et al. (2014), Koposov et al. (2015), and Li et al. (2016). Earlier studies hesitated to assign this structure to the Sgr stream, because simulations did not show any counterparts to these distant stars (e.g., Ruhland et al. 2011; Drake et al. 2013). Importantly, our model correctly replicates both the measured apocentric distance of 102.5 ± 2.5 kpc and the position angle found by Belokurov et al. (2014). Two main mechanisms have been invoked in the literature to explain the difference in apocentric distances of ∼50 kpc measured by Belokurov et al. (2014). According to Chakrabarti et al. (2014), the disparity is due to dynamical friction modulating the eccentricity of the orbit of Sgr. Given the regime of high initial separation and high progenitor mass investigated in our simulation, dynamical friction plays an important role in shrinking the pericentric and apocentric distances reached by the Sgr dwarf. On the other hand, Gibbons et al. (2014) find that the stream behavior depends primarily on the host potential and secondarily on the properties of the satellite progenitor.

Figure 10 shows the present-day distribution of the simulated Sgr stars in different projections of phase space. The panels on the left feature all Sgr star particles inside a Galactocentric radius of 100 kpc, along with stellar tracer data from 2MASS from Majewski et al. (2004). Overall, the simulation produces a reasonable match to the observed tidal debris characteristics shown in black squares. Unlike the stellar distribution presented in Gómez et al. (2015), we successfully reproduce the distant branch of the leading arm observed at R.A. ∼ 200° and distances >50 kpc. The model, however, predicts comparatively fewer stars at smaller distances. This comparison is limited, however, by the low number of Sgr stellar particles in the simulation and the resulting relatively sparse sampling. Similarly, the simulated data at R.A. ∼ 150° and decl. ∼ 20° are too coarse to distinguish the presence of the bifurcation in the leading stream detected by Belokurov et al. (2006). Observational and theoretical studies (e.g., Peñarrubia et al. 2010, 2011; Belokurov et al. 2014) have not yet reached a consensus regarding the dependence of the bifurcated tails on the morphology of the MW and the Sgr progenitor. More comprehensive simulations and detailed testing of galactic model parameters are needed in order to shed light on the origins of the stream bifurcation.

Zoom In Zoom Out Reset image size

The simulated stellar locations in (R.A., decl.) feature a systematic offset compared to the data, and the location of highest particle density does not line up exactly with the Sgr centroid. Such coherent shifts of the simulated stream occasionally appeared across the numerous test simulations carried out for this study. The combination of coarse sampling of possible Sun locations on the Fibonacci sphere (see Section 2.2) as well as the time interval of 25 Myr between subsequent output snapshots may be responsible for these offsets. Finer sampling and tuning of the simulation parameters could perhaps remove this issue. Alternatively, this angular offset may indicate that the simple spherical potential used to model the MW halo in this simulation is insufficient.

Additionally, the run of line-of-sight velocities for the leading arm (250° ≲ R.A. ≲ 150°) features a discrepancy with the data similar to that discussed in other studies (e.g., Helmi 2004; Law et al. 2005). Helmi (2004) suggests that a prolate model of the MW potential is required to solve this mismatch. However, Law et al. (2005) and Johnston et al. (2005) point out that such models introduce a discordance with the precession angles measured from stream M-giants. While these studies disagree on the exact nature of the MW triaxiality (oblate versus prolate), there is a consensus that it is difficult to reproduce the radial velocities of the leading arm with a simple spherically symmetric model (such as the one used in our study). On the other hand, some papers (Law & Majewski 2010; Gómez et al. 2015) have shown that the LMC can introduce significant perturbations to the phase-space distribution of Sgr debris. These hypotheses could be tested in future iterations of the orbital model presented in our study.

The set of panels on the right of Figure 10 show the predicted distribution of stars beyond distances of 100 kpc, where no Sgr stream debris have yet been identified with certainty. The large-scale tidal features labeled in Figure 8 appear as prominent clouds of stars in these panels. The structures located at R.A. ∼ 100°, ∼5°, and ∼250° correspond to those labeled N_W, N_E, and S in Figure 8, respectively. These distant branches can be distinguished from the closer debris in that area of the sky through line-of-sight velocities (color-coded with a different range for the distant stars), which are generally lower by ∼50 km s⁻¹ because they approach turnaround.

While no tidal streams have so far been identified at such large distances (Drake et al. 2013), Figure 10 includes the parameters of the few MW stars found at distances >100 kpc. A variety of stellar tracers have been used to map the inner halo. However, due to faint limits of the order of r ≲ 21, surveys so far have yielded only a dozen or so stars beyond 100 kpc. Deason et al. (2012) provide a sample of distant blue horizontal-branch and N-type carbon stars, later complemented by the M-giants detected by Bochanski et al. (2014). In their data from the UKIRT Infrared Deep Sky Survey (UKIDSS), Bochanski et al. (2014) were able to identify two stars with distances above 200 kpc, making them the most distant known MW stars. Figure 10 shows that these most distant stars appear consistent with the remote northern stream branch predicted in our simulation. Bochanski et al. (2014) already associated the M-giants in a distance range of 20–90 kpc in the UKIDSS data with Sgr. However, the more distant detections in their sample were not connected with the stream because the models of Law & Majewski (2010) do not feature stars beyond ∼75 kpc. Our study shows that some of the most distant known MW stars could have originated in Sgr.