THE SAGITTARIUS DWARF GALAXY: A MODEL FOR EVOLUTION IN A TRIAXIAL MILKY WAY HALO

We adopt a basic formalism similar to that described by LJM05 and LMJ09. The Milky Way is described by a smooth fixed gravitational potential consisting of a Miyamoto & Nagai (1975) disk, Hernquist spheroid, and a logarithmic halo. The gravitational potential is therefore specified directly through the isovelocity ellipsoids both for consistency with previous studies and because these are more analytically accessible than working with an underlying mass distribution from which the gravitational acceleration must be numerically derived.

The respective contribution of these components is given by

$$\begin{equation} \Phi _{\rm disk}=- \alpha {GM_{\rm disk} \over \sqrt{R^{2}+(a+\sqrt{z^{2}+b^{2}})^{2}}}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \Phi _{\rm sphere}=-{GM_{\rm sphere} \over r+c}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \Phi _{\rm halo}=v_{\rm halo}^2 \ln \big(C_1 x^2 + C_2 y^2 +C_3 x y + (z/q_z)^2 + r_{\rm halo}^2\big), \end{equation} \tag{ 3 }$$

where R/r are cylindrical/spherical radii respectively, and the various constants C₁, C₂, andC₃ are given by

$$\begin{equation} C_1 = \left(\frac{\hbox{cos}^2 \phi }{q_{1}^2} + \frac{\hbox{sin}^2 \phi }{q_{2}^2}\right), \end{equation} \tag{ 4 }$$

$$\begin{equation} C_2 = \left(\frac{\hbox{cos}^2 \phi }{q_{2}^2} + \frac{\hbox{sin}^2 \phi }{q_{1}^2}\right), \end{equation} \tag{ 5 }$$

$$\begin{equation} C_3 = 2 \, \hbox{sin}\phi \, \hbox{cos} \phi \left(\frac{1}{q_{1}^2} - \frac{1}{q_{2}^2}\right). \end{equation} \tag{ 6 }$$

As discussed by LMJ09, this form for the dark halo potential describes an ellipsoid rotated by an angle ϕ about the Galactic Z-axis, in which q₁ and q₂ are the axial flattenings along the equatorial axes and q_z is the axial flattening perpendicular to the Galactic disk. We adopt the convention that ϕ = 0° corresponds to q₁/q₂ coincident with the Galactic X/Y-axes respectively, and increases in the direction of positive Galactic longitude. Since only the ratios between the various q are important (as increasing/decreasing q₁/q₂/q_z in tandem is largely degenerate with the total halo mass normalization v_halo), we fix q₂ = 1.0.

We held fixed any parameter that was found not to be significantly constrained by Sgr (i.e., that did not significantly change the path of Sgr debris when varied). Following Johnston et al. (1999) and LJM05, we take M_disk = 1.0 × 10¹¹ M_☉, M_sphere = 3.4 × 10¹⁰ M_☉, a = 6.5 kpc, b = 0.26 kpc, and c = 0.7 kpc. In LJM05, we explored the possibility of scaling the total mass of the Galactic disk through the parameter α and found no compelling reason to adopt a value other than α = 1.0. We fix the distance to the Galactic center to be R_☉ = 8 kpc (e.g., Reid 1993; Groenewegen et al. 2008); in LJM05, we explored the possibility of varying R_☉, and found that for cases where the Sun lies approximately in the plane of the orbit (as is the case for Sgr) any reasonable choice of R_☉ does not have a particularly noticeable effect on the quality of fit of our simulations. Additionally, Figure 4 of LJM05 suggests that the ideal choice of the scale length of the Galactic halo (r_halo) is reasonably insensitive to the axial flattening adopted, and that a single fixed value of r_halo = 12 kpc should suffice. The normalization of the dark halo mass via the scale parameter v_halo is therefore entirely specified (for a given choice of q₁, q₂, and ϕ) by the requirement that the speed of the local standard of rest be v_LSR = 220 km s⁻¹ . As for R_☉, v_LSR is at present relatively unconstrained by observations of the Sgr stream (LJM05; although see discussion by Majewski et al. 2006; J. L. Carlin et al. 2010, in preparation) and we simply fix v_LSR to this conventional value. We note that if we were instead to adopt the proper motion of Sgr A* measured by Reid & Brunthaler (2004) it would imply v_LSR = 235.6 km s⁻¹ for our adopted value of R_☉ = 8.0 kpc.

In a recent contribution (LMJ09) we integrated massless test particles representing the Sgr dwarf along orbits in the Galactic potential defined by Equations (1)–(3) for a range⁴ of parameters ϕ = 0°–180°, q₁ = 0.6–1.8, and q_z = 0.6–1.8. While the accuracy of such simple orbit models is necessarily limited because actual tidal debris from a massive satellite will deviate slightly from these orbits, LMJ09 were able to identify a region of parameter space around ϕ ≈ 90°, q₁ ≈ 1.5, and q_z ≈ 1.25 generally favorable to producing a model satellite capable of matching the observational constraints enumerated in Section 2.2. We note that while such test-particle orbits technically identify a second similarly acceptable region of parameter space around ϕ ≈ 0°, q₁ ≈ 0.67, and q_z ≈ 0.83, the physical axis ratios of this region are degenerate with the first: Both have the short axis approximately aligned with the Galactic X-axis, and the longest with the Galactic Y. We choose to focus on the q_z>1 branch of the solution to avoid introducing degeneracy into our simulations.

Motivated by the successes of LMJ09, we perform N-body simulations in Galactic models with a range of parameters ϕ = 75°–115°, q₁ = 1.1–1.6, and q_z = 1.1–1.5. In each case, Sgr is integrated ∼8 Gyr backward in time along a simple test-particle orbit, and a 10⁵-particle model for the dwarf is integrated forward again in time from this point to the present day. The particles in our model of Sgr are initially distributed to generate a Plummer (1911) model

$$\begin{equation} \Phi _{\rm Sgr}=-{GM_{\rm Sgr, 0} \over \sqrt{r^2+r_{\rm 0}^2}}, \end{equation} \tag{ 7 }$$

where M_Sgr,0 is the initial mass of Sgr and r₀ = (M_Sgr,0/10⁹ M_☉)^1/3 kpc is its scale length. We have assumed spherical symmetry and isotropy of the orbits within the initial Sgr dwarf since: (1) there is no evidence yet of significant rotation in the main body of the satellite and (2) debris resulting from isotropic orbits will cover all phase-space regions accessible to non-isotropic satellites. The simulation particles are not associated with specific dark or light components, but instead trace the total mass of the satellite. The mutual interactions between these self-gravitating particles⁵ in the external Milky Way potential are calculated using a self-consistent field code (Hernquist & Ostriker 1992).

We summarize below the observational constraints that we adopt on the physical properties of the Milky Way–Sgr system, each of which is discussed in detail in Sections 3 and 4. Generally, we have adopted these constraints because we consider them to be the strongest, most unambiguous available. Note that we opt not to constrain the model to match the apparent bifurcation in the Sgr leading arm described by Belokurov et al. (2006) since recent observational data presented by Yanny et al. (2009) suggest that feature may arise due to substructure within the Sgr dwarf rather than the shape of the Galactic gravitational potential (see discussion in Sections 4.2 and 5.2).

While individual constraints are discussed in detail in Section 3, we summarize the following below.

• 1.  
  Sgr must lie at the angular coordinates, distance, and radial velocity described in Section 3.1.
• 2.  
  Sgr must have a velocity vector that lies within the orbital plane defined by recent trailing tidal debris.
• 3.  
  The radial velocity dispersion among trailing tidal debris must match that observed by Monaco et al. (2007) using high-resolution echelle spectroscopy.
• 4.  
  The run of radial velocities with orbital longitude along trailing tidal debris must match data presented by Majewski et al. (2004).
• 5.  
  The run of radial velocities with orbital longitude along leading tidal debris must match data presented by Law et al. (2004; see also LJM05).
• 6.  
  The angular location of leading tidal debris must match the high surface-brightness stream observed by Belokurov et al. (2006) in the SDSS.
• 7.  
  The angular width of leading tidal debris must match that observed by Belokurov et al. (2006).

It turns out to be easiest to constrain the properties of the Sgr dwarf itself first (Section 3), since many of them are largely independent of the exact form of the Galactic potential (to within the range explored in this paper), and the rest may be constrained contingent upon a specific potential. We then use these simulations to discriminate between different realizations of the underlying Galactic potential in Section 4. We note that the only distance information explicitly incorporated in our fitting method is that used in setting the distance to the Sgr dwarf itself. Our analysis therefore relies only upon the most well-constrained observations of the tidal debris (namely, the angular coordinates and radial velocities) and is insensitive to photometric distance calibration uncertainties which have hindered confident interpretation of previous studies. While we do not fit for the distances to Sgr tidal debris explicitly, we shall see in Sections 5 and  7.3 that the final model is nonetheless an excellent match to the distances derived for the observational data using the most reasonable assumptions about distance calibrations.

We assume that the Sgr dwarf is located at Galactic coordinates (l, b) = (56, − 142) (Majewski et al. 2003). Estimates of the distance to Sgr typically range from 24to28 kpc (see Table 2 of Kunder & Chaboyer 2009 for a summary); we adopt the distance scale D_Sgr = 28 kpc from Siegel et al. (2007) based on the Hubble Space Telescope (HST)/Advanced Camera for Survey (ACS) analysis of the stellar populations in the Sgr core. In LJM05 (which adopted D_Sgr = 24 kpc), we explored the effect of varying D_Sgr and found it to be largely degenerate with other parameters (such as R_☉) and not particularly well constrained, although models generally favored higher values of D_Sgr/R_☉. We therefore simply fix D_Sgr = 28 kpc to minimize the number of free parameters in our model. This places the Sgr dwarf at Galactic coordinates⁶(X, Y, Z) = (19.0, 2.7, − 6.9) kpc. The radial velocity of Sgr is taken to be v_rad = 171 km s⁻¹ (Ibata et al. 1997), and the dwarf is assumed to be presently moving toward the Galactic plane (Irwin et al. 1996; Dinescu et al. 2005).

With the advantages of an internally homogeneous selection function, M-giants drawn from the 2MASS survey clearly define the orbital plane of the Sgr stream, especially those high-latitude M-giants free of thick disk contamination and spanning ∼150° across the south Galactic hemisphere (Majewski et al. 2003). This trailing stream is visible as it wraps from the Sgr core (Λ_☉ = 0°) through the south Galactic hemisphere toward the Galactic anticenter (Λ_☉ ∼ 150°)—near the position of the previous apocenter of the Sgr orbit (Λ_☉ ∼ 170°). The stream appears to dwindle rapidly at Λ_☉ ≳ 140°, however, largely because the increasing age and decreasing metallicity (e.g., Chou et al. 2007, 2010) of this older tidal debris results in a decline in the prevalence of the M-giant tracer population.

The trend of radial velocities along the trailing stream is well defined (Majewski et al. 2004), enabling us to use kinematical information to prune Galactic interlopers and select a sample of extremely high-likelihood Sgr trailing arm stars. We draw such a sample from the Majewski et al. (2004) catalog, including all stars in the range Λ_☉ = 0°–140° for which radial velocity data has been obtained.⁷ We fit a fourth-order polynomial to the radial velocity data as a function of orbital longitude, using an iterative 2.5σ rejection algorithm to produce our final clean sample of trailing arm stars (Figure 1, left panel).⁸ Note that we have not applied an explicit distance selection cut to the sample of stars plotted in Figure 1; if we were to do so it would remove a majority of the crosses from the left-hand panel as these stars lie significantly closer than the Sgr trailing arm stream (see, e.g., Figure 2 of Majewski et al. 2004).

Zoom In Zoom Out Reset image size

The visible and well-defined portion of the trailing arm seen in the southern Galactic hemisphere 2MASS M-giants is particularly important because it is dynamically young (torn off in the last ∼0–2 pericentric passages, depending on the assumed mass of Sgr) and has therefore experienced very little differential precession from the main body of Sgr (e.g., Helmi 2004; Johnston et al. 2005). Observations (Majewski et al. 2003) indicate that the trailing debris stream is well-confined to a single plane, and numerical simulations show that N-body debris in the range Λ_☉ = 0°–140° have angular momentum vectors aligned with that of the present satellite to within 1 degree for a wide range of models of the Galactic gravitational potential (see, e.g., Figures 4 and 5 of Johnston et al. 2005).

The lack of orbital precession in the trailing stream allows us to determine robustly the instantaneous angular momentum vector of Sgr. Using the sample of trailing arm stars defined above, we use a GC3 cell count technique (Johnston et al. 1996) to define the plane that best fits all of the trailing arm stars and the present position of the Sgr dwarf. For any given orbital pole (l_p, b_p) the distance D_i of star i from the plane defined by the pole is

$$\begin{equation} D_i = \hat{n} \cdot (\vec{r}_i - \vec{r}_{\rm Sgr}), \end{equation} \tag{ 8 }$$

where $\hat{n}=(\textnormal {cos}\,b_{\rm p} \, \textnormal {cos} \,l_{\rm p}, \textnormal {cos}\,b_{\rm p} \, \textnormal {sin}\,l_{\rm p}, \textnormal {sin}\,b_{\rm p})$ is the orbital pole vector, $\vec{r}_i$ is the vector from the Galactic center to star i, and $\vec{r}_{\rm Sgr}$ is the vector from the Galactic center to Sgr.

The best-fitting plane is that which minimizes the sum D² = Σ_iD²_i over the 108 data points in the clean sample of trailing arm stars shown in Figure 1. In Figure 2, we plot D as a function of (l_p, b_p), and demonstrate that the best-fitting orbital pole is (l_p, b_p) = (2738, − 145). Unsurprisingly (since both the sample of stars and the technique were very similar), this is nearly identical to the trailing arm pole calculated by Johnston et al. (2005), although it differs slightly from the pole adopted by LJM05, who used the pole derived by Majewski et al. (2003) from fitting to both leading and trailing Sgr tidal debris.

Zoom In Zoom Out Reset image size

Since the radial velocity vector $\vec{v}_{\rm rad}$ and orbital pole (Section 3.2) of Sgr are known, the three-dimensional space velocity of Sgr ($\vec{v}_{\rm Sgr} = \vec{v}_{\rm rad} + \vec{v}_{\rm tan}$, where by definition $\vec{v}_{\rm rad} \cdot \vec{v}_{\rm tan} = 0$) can be uniquely specified by the magnitude v_tan of the velocity of Sgr tangential to the line of sight:

$$\begin{equation} \vec{v}_{\rm Sgr} \cdot \hat{n} = \vec{v}_{\rm rad} \cdot \hat{n} + \vec{v}_{\rm tan} \cdot \hat{n} = 0. \end{equation} \tag{ 9 }$$

In particular, for each choice of the Galactic halo parameters ϕ, q₁, q_z there is a value of v_tan that produces a trailing tidal stream whose velocities best match those shown in Figure 1. We converge upon this value in a manner similar to LJM05 by minimizing the χ² statistic:

$$\begin{equation} \chi _{v,{\rm trail}}^2 = \frac{1}{n_{v,{\rm trail}} - 3} \sum _i \frac{(v_{\rm model} [i] - v_{\rm obs,trail} [i])^2}{\sigma _{v}^2}, \end{equation} \tag{ 10 }$$

where v_model[i] is determined via linear interpolation of the model to the Λ_☉[i] of each of the n_v,trail = 108 velocity data points v_obs,trail[i].

There are two components governing the apparent velocity width of the stream; physical dispersion and observational uncertainty in the derived radial velocities. Majewski et al. (2004) find a net dispersion of 12 km s⁻¹ , composed of 6 km s⁻¹ intrinsic uncertainty in the velocity measurements and an underlying 10 km s⁻¹ physical dispersion in the stream. For the purposes of Equation (10), it is irrelevant what the source of the dispersion is (since both will have identical effects governing the offset of stars from the orbital path) and we take σ_v = 12 km s⁻¹ . Rather than perform computationally intensive N-body simulations to constrain v_tan for every choice of ϕ, q₁, q_z, we instead use test-particle orbits (see, e.g., LMJ09) to estimate the necessary value. Based on ∼150 N-body simulations with a wide range of Sgr masses and forms for the Galactic potential, we find a typical correction factor where the best-fit v_tan for N-body models is systematically 4 km s⁻¹ slower than that implied by χ² minimization using such test-particle orbits.

Considering the example case where ϕ = 90° and q₁ = 1.0, in Figure 3 (bottom panel) we illustrate the dependence of χ²_v,trail on v_tan for three different values of q_z. Note that there is a well-defined best value (which we call v_best) of v_tan for each q_z; we estimate the location of the minimum by fitting a simple third-order polynomial to χ²_v,trail as a function of v_tan. The resulting values of v_best are well behaved as a function of q_z, as illustrated in the top panel of Figure 3.

Zoom In Zoom Out Reset image size

For each choice of ϕ and q₁, it is possible to construct a similar relation. For ϕ = 90°, the relation between v_best and q_z is practically unchanged (to within ∼1–2 km s⁻¹ ) for any choice of q₁ = 0.8–1.6. For other choices of ϕ and q₁, the deviation of the relation from that shown in Figure 3 grows as ϕ departs from 90°, but for the range of ϕ = 75°–115° considered here these relations differ by less than 10 km s⁻¹ .

Data in the longitude range Λ_☉ = 25°–90° consist of a relatively orderly sample of trailing arm stars as they pass through the south Galactic cap moving predominantly transverse to our line of sight, and stars in this range can provide the most accurate measurement of the intrinsic velocity dispersion of the stream.⁹ Majewski et al. (2004) performed such an analysis, concluding that for a reasonable estimate of the observational uncertainty that the intrinsic velocity dispersion of the stream is σ_v,obs = 10.0 ± 1.7 km s⁻¹ . More recent, higher spectral resolution observations¹⁰ by Monaco et al. (2007) have revised this value slightly to σ_v,obs = 8.3 ± 0.9 km s⁻¹ .

Since the velocity dispersion in the tidal debris streams reflects the mass present within the tidal radius of Sgr on the orbit immediately prior to that debris becoming unbound (see, e.g., Johnston et al. 1996; Johnston 1998; Peñarrubia et al. 2008), it can therefore be used to constrain the mass of the dSph (e.g., LJM05). We performed a series of ∼160 N-body simulations in which satellites with 10 different initial masses ranging from M_Sgr,0 = 6.67 × 10⁷to 2.57 × 10⁹ M_☉ were integrated along reasonable orbits for various lengths of time in different realizations of the Milky Way gravitational potential. The resulting relation between the final bound Sgr mass M_Sgr and the velocity dispersion σ_v is shown in Figure 4. Extrapolating from the best-fit polynomial to the curve, we conclude that a present bound mass M_Sgr = 2.5^+1.3_−1.0 × 10⁸ M_☉ produces the best fit to the observed velocity dispersion. Adopting a typical mass-loss history (generally similar along reasonable orbits in the range of Galactic potentials explored here), this implies that the original mass of Sgr was M_Sgr,0 = 6.4^+3.6_−2.4 × 10⁸ M_☉ if it has been orbiting in the Galactic potential for ∼8 Gyr. We adopt M_Sgr,0 = 6.4 × 10⁸ M_☉, corresponding to an initial Sgr scale length r₀ = 0.85 kpc (see Equation (7)).

Zoom In Zoom Out Reset image size

As defined in Section 3.3, we use the statistic χ_v,trail to quantify the fit of a given N-body simulation to the trailing arm velocity trend. While χ_v,trail has already been minimized for a given choice of ϕ, q₁, q_z by determining the orbital velocity of the satellite, the precise minimum value can vary slightly between different halo realizations and we therefore also use it as a constraint on the properties of the Milky Way.

Although both the leading and trailing tidal tails are traced by the M-giant population, the location of Sgr slightly past pericenter along its orbit in combination with foreshortening of the orbital rosettes caused by our location 8 kpc from the center of the Milky Way means that while the trailing arm extends for ∼150° across the southern Galactic hemisphere the leading arm is clearly visible only for ∼75° (a significant section of which is obscured as it passes behind the Galactic disk/bulge). The exact track of leading Sgr debris as it passes through the north Galactic cap toward the anticenter is therefore uncertain in the 2MASS data since the M-giant tracer population is no longer as prominent in this older section of the stream. In contrast, SDSS imaging is sensitive to older Sgr populations, and clearly shows (Belokurov et al. 2006; Newberg et al. 2007; Cole et al. 2008) the Sgr leading arm arcing through the northern Galactic cap toward re-entry into the Galactic disk in the direction of the anticenter.

Belokurov et al. (2006) initially suggested that there may be a bifurcation in the angular position of the Sgr stream in the northern Galactic hemisphere, with a high surface-brightness "A" stream representing the primary wrap of the leading tidal arm passing through the coordinates (α, δ) ∼ (160°, 20°) and a lower surface-brightness "B" stream representing a secondary wrap of the trailing tidal arm offset to higher declination around (α, δ) ∼ (160°, 30°) (see also Fellhauer et al. 2006). We adopt the positions of the "A"-stream survey fields of Belokurov et al. (2006) along this high surface-brightness branch as our constraint on the orbital path of leading Sgr tidal debris. These 17 pointings are tabulated in Table 1. Note that Belokurov et al. (2006) adopted a photometric distance scale calibrated against a distance to Sgr of 25 kpc. The distances listed in Table 1 have been rescaled to our adopted D_Sgr = 28 kpc by multiplying by a factor of 1.12.

We note that the dynamical origin of the lower surface-brightness B stream is uncertain: recent observational data presented by Yanny et al. (2009) indicate that the A and B streams have indistinguishable trends in distance, radial velocity, metallicity, and mix of stellar types. This suggests that the A and B branches represent debris stripped from Sgr at similar times, and that the bifurcation may likely be within a single phase of tidal debris, originating in substructure within the initial Sgr dwarf (see discussion in Section 5.2). Since we do not attempt to constrain the internal structure of the Sgr dwarf (we focus here on exploring the effects of triaxiality in the gravitational potential of the Milky Way, and such an effort is beyond the scope of this contribution), we therefore do not explicitly constrain our model to match the B branch of the bifurcation.

Numerically, we define the statistic

$$\begin{equation} \chi _{\delta }^2 = \frac{1}{n_{\delta } - 3} \sum _i \frac{(\delta _{\rm model}[i] - \delta _{\rm obs}[i])^2}{\sigma _{\delta }^2} \end{equation} \tag{ 11 }$$

describing the accuracy with which the declination of the N-body model traces the A stream of Belokurov et al. (2006), and

$$\begin{equation} \chi _{w}^2 = \frac{(w_{\delta, {\rm model}} - w_{\delta, {\rm obs}})^2}{\sigma _{w}^2} \end{equation} \tag{ 12 }$$

describing the accuracy with which the angular width of the N-body model matches the SDSS data at a right ascension α = 180° (the non-spherical gravitational potential gives rise to angular broadening of the leading debris stream in this region due to differential precession). Based on the data presented by Belokurov et al. (2006), we estimate that the 1σ angular width of the stream as seen in SDSS is w_δ,obs = 52, with an uncertainty σ_w = 13. In contrast to Equation (10) (which used Λ_☉[i] as the variable coordinate), δ_model[i] is determined via linear interpolation of the model to the right ascension α[i] of each of the n_δ = 17 SDSS survey fields along the main branch of the Sgr stream, each of which is assumed to have an angular uncertainty,¹¹σ_δ = 19.

While the trend of leading tidal debris is difficult to trace with 2MASS imaging data alone, the stream can be recovered using spectroscopy to identify the radial velocity signature of the stream. In Figure 1, we show our selection function for stars in the Sgr leading arm with K_s>9 (see Section 7.3) from the observational data presented by Law et al. (2004; see also LJM05), using a method similar to that adopted for trailing arm stars. This velocity trend has recently been confirmed by spectroscopic observation of a larger number of stars drawn from the SDSS imaging of the leading tidal stream (Yanny et al. 2009), in contrast with the preliminary conjecture by Newberg et al. (2007) that the velocity trend may be artificially produced by confusion in the leading arm M-giant sample. This gives our final constraint on the stream

$$\begin{equation} \chi _{v,{\rm lead}}^2 = \frac{1}{n_{v,{\rm lead}} - 3} \sum _i \frac{(v_{\rm model} [i] - v_{\rm obs,lead} [i])^2}{\sigma _{v}^2}, \end{equation} \tag{ 13 }$$

where v_model[i] is determined via linear interpolation of the model to the orbital longitude Λ_☉[i] of each of the n_v,lead = 94 observational data points in the 2MASS leading arm radial velocity sample. As for the trailing arm velocity constraint, we adopt an uncertainty σ_v = 12 km s⁻¹ for each radial velocity measurement representing a combination of observational uncertainty and intrinsic stream width.

In general, the parameters ϕ, q₁, q_z are not separable and must be constrained simultaneously by χ_v,trail, χ_v,lead, χ_δ, and χ_w. We combine these four constraints into a single statistic:

$$\begin{equation} \chi ^2 = \chi _{\rm v, trail}^2 + \chi _{\rm v, lead}^2 + \chi _{\delta }^2 + \chi _w^2 \end{equation} \tag{ 14 }$$

and seek the combination of ϕ, q₁, q_z that gives the global minimum for this combined statistic.

We performed a series of ∼500 N-body simulations covering the range ϕ = 75°–115°, q₁ = 1.1–1.6, q_z = 1.1–1.5 with a grid spacing of 5° in ϕ, 0.05 in q_z, and 0.1 in q₁. The resulting values of χ as a function of ϕ, q₁, q_z are shown in Figure 5 and demonstrate a general preference for values of q₁ ≈ q_z in the range ϕ ≈ 90°–105°. Interpolating along the grid surface, the optimal fit (χ = 3.4) is achieved for ϕ = 97°, q₁ = 1.38, and q_z = 1.36. We use these parameters for our final best-fit N-body simulation.

Zoom In Zoom Out Reset image size

Since this optimal fit has χ>1 there are obviously some residual imperfections in the model because different observational constraints favor different regions of parameter space. The greatest mismatch lies in the leading arm radial velocities (χ_v,lead = 2.7), followed by the leading arm angular coordinates (χ_δ = 1.6), the trailing arm radial velocities (χ_v,trail = 1.3), and lastly the angular width of the leading stream (χ_w = 0.9). For comparison, the spherical halo model discussed by LJM05 is clearly a poor fit to the observational data (see Figure 10 of LJM05) and has χ ≈ 9. In general, we deem simulations with χ ≳ 5 to be visibly poor fits.

The optimal shape for the Galactic dark halo therefore appears to be an approximately oblate ellipsoid whose short axis is located in the Galactic disk along the direction (l, b) = (7°, 0°). We note that this differs somewhat from that derived by LMJ09, who found a more strongly triaxial ellipsoid with ϕ = 90°, q₁ = 1.50, and q_z = 1.25 by matching observed tidal debris to the orbital path of a point-mass orbiting in the Galactic potential. It is well known, however, that actual leading/trailing arm debris falls inside/outside of the orbit, respectively (e.g., Eyre & Binney 2009), meaning that a slightly different shape will be derived for the Galactic halo when attempting to match simulated debris to observed debris (as done in this contribution), rather than a simulated orbit to observed debris (as in LMJ09).

In LJM05, we presented possible models for the tidal disruption of Sgr under three different assumptions; that the Galactic halo was axisymmetric in the disk plane and either oblate, spherical, or prolate in the direction perpendicular to the disk. Neither of these three models were entirely satisfactory as they could not simultaneously reproduce the angular position, distance, and radial velocity trends of the leading arm. In contrast, our new model,¹² orbiting in a non-axisymmetric dark matter halo, reproduces all of these constraints to reasonable accuracy as demonstrated by Figure 6. Indeed, we note that even though we did not explicitly fit to the distances of the SDSS leading arm data from Belokurov et al. (2006), the new model nonetheless matches them to within the uncertainty of the photometric calibration. While it is not possible to test the agreement of this new model with the distances of the M-giant sample directly (due to the evolution of the MDF along the Sgr tidal streams), we demonstrate in Section 7.3 that the model satisfactorily reproduces their apparent magnitude trend.

Zoom In Zoom Out Reset image size

Since the Sgr streams frequently overlap themselves in all of the observable coordinates, it is useful to isolate specific wraps of the leading/trailing streams in order to discuss individual wraps more easily. By sorting the streams according to angular momentum (relative to the Sgr core) and time at which a given debris particle became unbound from the model satellite this is easily accomplished in the N-body simulations. We color-code debris particles according to the time they became unbound from the satellite: black points are currently bound to Sgr, orange points were stripped from the dwarf during the last two perigalactic passages (0–1.3 Gyr ago), magenta points during the previous two (1.3–3.2 Gyr ago), etc. (see Figure 7 for a detailed accounting of the color-coding scheme). As a convenient shorthand notation, we describe specific wraps of the streams as follows (see summary in Figure 8).

• L1. Primary wrap of the leading arm of Sgr (i.e., 0°–360° angular separation from the dwarf). Traced by orange/magenta points near the Sgr dwarf, magenta/cyan/green points at larger angular separations.
• L2. Secondary wrap of the leading arm of Sgr (i.e., 360°–720° angular separation from the dwarf). Traced almost entirely by cyan/green points.
• T1. Primary wrap of the trailing arm of Sgr (i.e., 0°–360° angular separation from the dwarf). Traced by orange/magenta points near the Sgr dwarf, magenta/cyan/green points at larger angular separations.
• T1. Secondary wrap of the trailing arm of Sgr (i.e., 360°–720° angular separation from the dwarf). Traced almost entirely by cyan/green points.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The major observational distinction of the present model from the oblate-halo model discussed by LJM05 is that the primary L1 wrap of the leading Sgr tidal stream is no longer predicted to pass closer than ∼10 kpc to the Sun; the center of the L1 streamer intersects the Galactic disk at a distance of 13.5 kpc in the direction of the Galactic anticenter (see, e.g., Figures 6 and 8). This is consistent with the absence of net vertical flow observed in the solar neighborhood in recent radial velocities surveys (Seabroke et al. 2008).

Due to the triaxial Galactic potential, the orbital period and pericenter/apocenter of Sgr varies from orbit to orbit (Figure 7); the mean values of these parameters are 0.93 Gyr and 19/59 kpc, respectively, over the past 8 Gyr. The present space velocity of Sgr in this model is (U, V, W) = (230, −35, 195) km s⁻¹ . For our assumptions about the distance to Sgr, the local standard of rest, and the solar peculiar motion this translates to a proper motion of (μ_lcosb, μ_b) = (−2.16, 1.73) mas yr⁻¹ ([μ_αcosδ, μ_δ] = [ − 2.45, − 1.30] mas yr⁻¹ in celestial coordinates) in the heliocentric rest frame. This agrees within 1.9σ with the proper motion measurements of Dinescu et al. (2005) who found (μ_lcosb, μ_b) = (−2.35 ±  0.20, 2.07 ±  0.20) mas yr⁻¹ based on Southern Proper Motion Catalog values for 2MASS-selected giant star candidates, and within 2.2σ of the values (μ_lcosb, μ_b) = (−2.62 ± 0.22, 1.87 ± 0.19) mas yr⁻¹ determined by Pryor et al. (2010). Our tangential velocity for Sgr (v_tan = 251 km s⁻¹ ) is also in good agreement with the measures of Irwin et al. (1996; v_tan = 250 ± 90 km s⁻¹ ) and Ibata et al. (2001a; v_tan = 280 ± 20 km s⁻¹ ). The mass of the Milky Way within 50 kpc in this model is 4.4 × 10¹¹ M_☉ (i.e., slightly smaller than that obtained by Sakamoto et al. 2003), corresponding to a halo scaling parameter v_halo = 121.9 km s⁻¹ (see Equation (3)).

The present mass of the Sgr dwarf in this model is M_Sgr = 2.5^+1.3_−1.0 × 10⁸ M_☉; taking the apparent magnitude of the bound core of Sgr to be m_V = 3.63 (Majewski et al. 2003) gives an absolute magnitude M_V = −13.64 and luminosity L_Sgr = 2.4 × 10⁷ L_☉ for our adopted Sgr distance modulus ([m − M]₀ = 17.27; Siegel et al. 2007) implying a mass/light ratio of M_Sgr/L_Sgr = 10⁺⁶₋₄ M_☉/L_☉. This is at the low end of the range M_Sgr/L_Sgr = 14 − 36 M_☉/L_☉ estimated by LJM05; this difference reflects both our greater adopted distance to Sgr and our choice to use the higher-resolution Monaco et al. (2007) measurement of the dispersion of the Sgr stream instead of the previous Majewski et al. (2004) value. We note, however, that the derived mass is somewhat dependent on the structural model adopted for the Sgr dwarf. At present, the total interaction time of Sgr (and hence the initial mass of the satellite) is relatively unconstrained: with reasonable assumptions for its initial mass the length of the tidal streams observed in 2MASS and SDSS suggest that Sgr must have been disrupting for at least the last 3 Gyr (roughly corresponding to orange/magenta-colored points in Figure 6). We have opted to run our simulation for a significantly longer time (∼8 Gyr) in order to make concrete predictions about where Sgr debris from earlier epochs may be expected; observational confirmation of the presence or absence of such older debris streams will provide a strong constraint on the actual length of the interaction time. Future models that more accurately treat the dark/light matter separation and are constrained to reproduce the bar-like core structure of Sgr (e.g., E. Łokas et al. 2010, in preparation) may provide more stringent constraints on the M/L ratio of the dwarf.

As indicated by Figure 6, our simulations match the high surface-brightness (hereafter HSB) track of the leading arm in the "Field of Streams" region described by Belokurov et al. (2006) but do not produce a bifurcation in the form of a low surface-brightness (hereafter LSB) branch offset to higher declination. Fellhauer et al. (2006) previously explained this bifurcation as the juxtaposition of the young (<4 Gyr) L1 wrap of the leading arm with the old (>7.4 Gyr) T2 wrap of the trailing arm.¹³ This explanation is not viable for the current model. As illustrated in Figure 9, the distance/radial velocity trends of the L1 (black points) and T2 (green points) streams are markedly different, in contrast to Yanny et al. (2009) who demonstrate that the distance and radial velocity trends of the LSB and HSB streams are indistinguishable.

Zoom In Zoom Out Reset image size

However, we note that the Yanny et al. (2009) observations are not reproduced well by the Fellhauer et al. (2006) model either: this model predicts a factor ∼2 difference in distance between the LSB and HSB streams at α = 120° that is not observed. Additionally, given the strong evolution in the MDF with increasing separation Λ_☉ from the Sgr dwarf (e.g., Chou et al. 2007) we should not expect the metallicity of the LSB and HSB streams to be similar (as observed by Yanny et al. 2009) if they represent extremely different dynamical ages of tidal debris. Certainly, if the LSB stream were significantly older than the HSB we should not expect the LSB stream to contain as many M-giants as the HSB stream since M-giants are a young tracer population.

N-body models of Sgr tidal destruction also predict that the young T1 trailing arm should pass through the Field of Streams region. In both the model presented here and the model described by Fellhauer et al. (2006), the T1 wrap traverses this region from α ∼ 120° to α ∼ 220° with a distance varying from ∼40 kpc to ∼20 kpc, and radial velocity varying from ∼125 km s⁻¹ to ∼−125 km s⁻¹ . Such tidal debris (gray points in Figure 9) is not observed in the Field of Streams (although c.f. Niederste-Ostholt et al. 2009 and discussion in D. R. Law et al. 2010, in preparation), suggesting that the Sgr tidal streams are insufficiently long to produce T1 trailing tidal debris in this region of sky (or at least, that the stellar density in the stream is declining sufficiently rapidly that it has not been detected in the SDSS). This places a limit on the length of the Sgr tidal stream similar to that derived from other observational constraints that the Sgr tails represent ≲3 Gyr worth of tidal debris. If the Sgr tails are so young, we should not expect any debris corresponding to earlier wraps (i.e., T2/L2) to exist and it is correspondingly impossible to interpret the bifurcation as the juxtaposition of Sgr tidal debris torn from Sgr at significantly different times.

We therefore favor the explanation proffered by Yanny et al. (2009) that both the LSB and HSB streams may represent debris stripped from Sgr at similar times. Indeed, we note that models of Sgr orbiting in a non-spherical Galactic potential exhibit precessional broadening of the angular width of the leading debris stream in this region which manifests as a LSB "spray" to higher declinations similar to that of the LSB stream (this effect is increasingly noticeable for halos with stronger departures from spherical symmetry). This presents an alternative scenario wherein the bifurcation may arise due to substructure within the initial Sgr dwarf, such as a smaller dwarf satellite that may have been bound to Sgr when it fell into the gravitational potential of the Milky Way (see, e.g., Chou et al. 2010; D. R. Law et al. 2010, in preparation), or simply to anisotropy of the stellar orbits within the original dwarf which might be expected to give rise to tidal debris occupying only a subset of the angular width covered by our isotropic N-body model.

The detailed nature of the LSB branch in the SDSS view of the Sgr stream is therefore unknown at present, and is not possible to reconcile with N-body models that fit the unambiguous primary leading and trailing tidal streams illustrated in Figure 6. While this feature may ultimately provide a constraint on the internal structure of the initial dSph satellite, exploration of such structure is beyond the scope of this contribution.

6.1.1. Tracing the Underlying Mass Distribution

In Section 2.1, we introduced flattening directly into the contours of the gravitational potential because this was analytically tractable and consistent with the method of previous studies (in particular LJM05). It is well known, however, that as the ellipticity of the potential becomes more extreme, the corresponding mass distribution required to generate it may become unphysical (i.e., have regions of negative mass). In order for any model of the Galactic halo to be believed, it must correspond to a physically realizable density distribution.

In Figure 10, we show the corresponding contours of the mass distribution underlying our best-fit model for the Milky Way (as calculated by taking the Laplacian of the gravitational potential). As expected, the mass distribution is less spherical than the isopotential surface; the contours of the dark matter halo are approximately elliptical at small radii (i.e., within ∼40 kpc) with (c/a)_ρ = 0.44, (b/a)_ρ = 0.97, and T_ρ = 0.07. Although some pinching of the dark matter density contours is apparent at larger radii the mass distribution is everywhere positive. Figure 10 also demonstrates that while addition of the baryonic mass component significantly changes the net density distribution at small radii and near the Galactic disk plane, the density distribution at larger radii probed by the orbit of Sgr is dominated by the dark halo component.

Zoom In Zoom Out Reset image size

We additionally confirmed the validity of our results by constructing a triaxial NFW (Navarro et al. 1996) halo model in which the axial flattenings were introduced in a more physically motivated manner directly into the mass distribution and where the resulting accelerations were computed using the approximate form of the corresponding potential given by Lee & Suto (2003). As expected, there are slight differences in the path of satellites orbiting in this NFW halo and in our best-fit logarithmic potential, but these differences are minor and likely due to the fact that we did not fine tune the NFW model.

6.1.2. Orientation of Principal Axes

In LJM05, we conducted an extensive parameter space search in an effort to find a combination of Milky Way + Sgr parameters capable of simultaneously matching both the observed angular position and distance/radial velocity trends of leading tidal debris. In particular, we explored the effect of varying the proper motion of Sgr (v_tan), the distance to Sgr (D_Sgr), the distance from the Sun to the Galactic center (R_☉), the total mass of the Galactic disk (via the scaling parameter α in Equation (1)), the dark halo scale length (d), the dark halo flattening perpendicular to the disk (q_z), and the total mass scale of the Milky Way as parameterized by the rotation speed of the local standard of rest (v_LSR).

Although the Galactic disk contributes significantly to the overall flattening of our adopted Galactic potential at small radii (e.g., Figure 10), it has only a minor effect at the range of radii (r ∼ 20–60 kpc) probed by the orbit of Sgr. As demonstrated by Figure 11, the dark halo in our model is the dominant contributor to the total gravitational potential at these radii. Indeed, if we decrease the mass of the Galactic disk by 50% (by setting α = 0.5 in Equation (1)) and increase the total halo mass correspondingly (so that the local standard of rest remains constant at v_LSR = 220 km s⁻¹ ), we can repeat the parameter space search described in Section 4. Scaling the proper motion v_tan of Sgr to again best fit the trailing arm velocity trend, we find optimal values for q₁, q_z, ϕ almost identical to those derived previously (although the overall χ² of the best-fit model is not as good as for the case α = 1.0). We are therefore confident that the details of the Galactic disk model do not significantly bias our results.

Zoom In Zoom Out Reset image size

Similarly, the details of the triaxial Galactic bar (which we did not include in our model) and stellar halo are unlikely to change our results. Neither are sufficiently massive to have a great effect on the gravitational potential at distances r ≳ 20 kpc, and neither have orientations akin to that which we derive for the dark halo. The major axis of the bar is thought to lie within ∼15°–20° of the X-axis (e.g., Blitz & Spergel 1991; Nakada et al. 1991; Morris & Serabyn 1996; Babusiaux & Gilmore 2005), similar to the minor axis of the dark halo. While the stellar halo is significantly triaxial ([c/a]_ρ ≈ 0.65, [b/a]_ρ = 0.75, T_ρ = 0.76; Newberg & Yanny 2006; see also Larsen & Humphreys 1996; Xu et al. 2006) with a major axis lying approximately in the Galactic plane (∼50°–70° from the X-axis) the minor axis is not apparently contained within the X–Y plane but located ∼13° from the Z-axis.

Throughout the entirety of the LJM05 parameter space search, the only parameter that was found to have a significant effect on the angular position/distance/radial velocity trend of the Sgr leading arm was the flattening of the dark halo q_z. Although no single choice of q_z can simultaneously reproduce all observational constraints, we have demonstrated (see also LJM09) that a fully triaxial model in which the short axis is approximately aligned with the Galactic X-axis is capable of doing so.

Such a triaxial halo is therefore obviously attractive in the sense that only N-body models in such a halo succeed in modeling the Sgr–Milky Way system where previous efforts (e.g., Helmi et al. 2004; LJM05; Fellhauer et al. 2006; Martínez-Delgado et al. 2007) have failed. Despite the extensive exploration of parameter space undertaken by ourselves (this paper; LJM05; LMJ09) and other groups (e.g., Helmi et al. 2004; Fellhauer et al. 2006; Martínez-Delgado et al. 2007) however, the space explored to date is far from exhaustive and it is possible that the triaxial halo model may prove to be simply a numerical crutch that mimics the effect of some as-yet unidentified alternative. In the following sections, we briefly discuss a few possibilities that may obviate the need for triaxiality to explain the extant data.

6.2.1. Alternative Formulations for the Gravitational Potential

6.2.2. Evolution in the Orbit

6.2.3. Dwarf–Tail Gravitational Interaction

6.2.4. Gravitational Influence of the Magellanic Clouds

The Large Magellanic Cloud (LMC) is ∼50 kpc distant in the direction (l, b) = (2805, − 328) (van der Marel et al. 2002), placing it within ∼0.5 kpc of the Galactic Y–Z plane. Based on great-circle fits to the Magellanic stream, the pole of the system appears to be (l_p, b_p) = (1885, − 75) (Nidever et al. 2008). That is, the pole of the Magellanic stream is aligned with the short axis of the Galactic dark halo as derived in this contribution to within ∼1° in Galactic longitude.¹⁶ This suggests two possibilities: (1) the Magellanic Clouds may have fallen into the Milky Way along the plane of the long and intermediate axes of the dark matter halo (i.e., the alignment is a consequence of the dark matter distribution) or (2) the Magellanic Clouds may exert a significant gravitational acceleration on the Sgr dwarf, with the mass of the Clouds at large radii in the Y–Z plane mimicking the effect of an extended dark matter distribution in this plane (i.e., the alignment is the cause of the apparent dark matter distribution).

A precise description of the second scenario would require a comprehensive model for the orbital and mass-loss history of the Magellanic Clouds, but it is possible to understand the magnitude of their gravitational influence on the path of Sgr tidal debris by integrating massless test particles along the Sgr orbit in a Milky Way with a spherical dark halo and simplified model of the LMC. We first test the effect of introducing a fixed point mass into the gravitational potential of the Milky Way at the present location of the LMC, whose mass is allowed to vary from 1%–10% of the total mass of the Milky Way within 50 kpc (i.e., M_LMC = 6 × 10⁹ M_☉ − 6 × 10¹⁰ M_☉). As illustrated in Figure 12, a fixed LMC with a mass 10% that of the Milky Way can produce a significant deviation on the order of 20° in the orbital path of leading Sgr tidal debris.

Zoom In Zoom Out Reset image size

Of course the LMC is not fixed in space, but describes an orbit about the Milky Way confined approximately to the Y–Z plane. We therefore develop a second, slightly more realistic model in which the LMC is described in a time-integrated sense by a ring of mass at a radius of 50 kpc in the Y–Z plane. The gravitational acceleration at an arbitrary point in space resulting from such a configuration is described using an analytical approximation to the numerically evaluated elliptical integrals (see derivation for the electrostatic case by Zypman 2006). This second model produces even more noticeable deviations from the original orbital path, with appreciable shifts in both the leading arm angular position and radial velocities expected for models in the range of masses 1%–10% that of the Milky Way.

The detailed effect of the Magellanic Clouds on the orbit of Sgr is extremely uncertain. Both our point-mass and ring-mass models are obviously oversimplifications, neglect the SMC entirely, and assume that the gravitational influence of the LMC has been fixed over the entire interaction history of the Sgr dwarf (∼8 Gyr in the current model). Since the Clouds are currently thought to be near the pericenter of their orbit, for the majority of the lifetime of Sgr they likely lay at greater distances where their gravitational force would be smaller. Indeed, some recent lines of investigation suggest that they may have only recently fallen into the gravitational potential of the Milky Way for the first time (e.g., Besla et al. 2007), drastically limiting the time during which they could have torqued the orbit of the Sgr dwarf. Of the various orbital paths shown in Figure 12, only the triaxial halo model actually succeeds in reproducing all of the observational data. However, given the myriad uncertainties in the Milky Way–LMC–Sgr system the important conclusion to draw is that gravitational perturbations from the Magellanic Clouds may be sufficient to torque the leading arm of the Sgr stream significantly, especially if the effect is enhanced by the response of the Milky Way dark matter halo itself to the passage of the Clouds (see, e.g., Weinberg 1998).

6.2.5. Alternative Gravity

Fortunately, the model presented here for the shape, size, and orientation of the Galactic gravitational potential is easily testable; it must be capable of explaining the orbital paths of other Galactic satellites that similarly provide constraints on the mass distribution in the Milky Way. Typically, previous studies of tidal debris systems in the Galactic halo have only explored the effects of changing the axial scale length perpendicular to the Galactic disk (i.e., q_z), and it is unclear whether satisfactory solutions for other satellites can be found in such a strongly non-axisymmetric potential as that suggested here. In particular, can realistic orbits be found to match such systems as the Monoceros tidal stream (e.g., Peñarrubia et al. 2005), the Magellanic Clouds¹⁷ (e.g., Besla et al. 2007; Rŭžička et al. 2007), the "orphan stream" (Belokurov et al. 2007), the GD-1 stream (Koposov et al. 2010), the newfound Cetus Polar Stream (Newberg et al. 2009), and tidally disrupting Galactic globular clusters such as Pal 5 (Odenkirchen et al. 2009)? Similarly, is a triaxial Galactic halo consistent with observations of the Galactic H i disk (e.g., Olling & Merrifield 2000; Kalberla et al. 2005; and references therein)? It is intriguing that Saha et al. (2009) discuss evidence to suggest that the dark matter halo must be similarly non-axisymmetric in the Galactic XY plane in order to explain the asymmetric flaring of the disk (although cf. Narayan et al. 2005). Future observations of hypervelocity stars at large distances may also provide key diagnostics (e.g., Gnedin et al. 2005; Yu & Madau 2007; Perets et al. 2009).

Generally, we expect that star formation in Sgr is likely to have occurred near the bottom of its gravitational potential well, and that stars thus formed in the core would gradually diffuse outward with time. Remnants of such stellar populations would tend to persist in the Sgr core until the present day and we therefore assume that the stellar populations present in the Sgr streams correspond to counterparts in the Sgr core. These Sgr core populations have been traced using high-precision HST-ACS photometry by Siegel et al. (2007), who find five major star formation episodes.

Population A. M54 metal-poor population.

[Fe/H] = −1.7 ± 0.1, age 13 ± 1 Gyr.

Population B. Sgr metal-poor population.

[Fe/H] = −1.2 ± 0.1, age 11 ± 1 Gyr.

Population C. Sgr intermediate population.

Three sub-populations:

C1: [Fe/H] = −0.65 ± 0.1, age 7 ± 0.6 Gyr,

C2: [Fe/H] = −0.45 ± 0.1, age 6 ± 0.6 Gyr,

C3: $\hbox{[Fe/H]}= -0.3\phantom{5} \pm 0.1$, age 5 ± 0.6 Gyr.

Population D. Sgr young population.

[Fe/H] = −0.05 ± 0.15, age 2.5 ± 0.5 Gyr.

Population E. Sgr very young population.

[Fe/H] = +0.55 ± 0.05, age 0.75 ± 0.1 Gyr.

Younger stellar populations might be realistically assumed to be more tightly bound than older populations, which have had more time to diffuse outward in the dwarf. Motivated by the presence of such population gradients within other Galactic dSphs (e.g., Tolstoy et al. 2004) for which the more metal-rich populations are more centrally concentrated and have lower velocity dispersion, we tag particles in our best-fit model of the Sgr dwarf to correspond to different stellar populations according to their total energy in the pre-interaction dSph. The energy of an individual particle in the satellite is given by

$$\begin{equation} E[i] = \frac{1}{2} v[i]^2 + \Phi [i], \end{equation} \tag{ 15 }$$

where v[i] is the velocity and Φ[i] the internal potential energy of particle i relative to the satellite. For convenience, we define E* to be an order-sorted version of the energy distribution, ranging from 0 (most strongly bound particle) to 1 (least bound particle) with a uniform distribution in the range 0–1.

We assign particles with a given range of E* to each population as illustrated in Figure 13, the locations of the energy divisions (and hence the total number of particles assigned to each burst of star formation) are free parameters. Once a particle has been assigned to a given population, we assign it a metallicity and age based on a Gaussian random selection about the mean and standard deviation of the population given in the list above, and randomly choose an infrared color uniformly from the range J − K_s = 1.0 − 1.1. Adopting Marigo et al. (2008) isochrones¹⁸ with a Chabrier (2001) IMF for each of the simulated populations, we calculate the absolute K_s magnitude for each simulated particle, adopting a 50%/50% mix of red giant branch (RGB) and asymptotic giant branch (AGB) stars. This particular mix of RGB/AGB stars is not intended to match any particular distribution, but to illustrate the general range of magnitudes expected for stars in each evolutionary phase. We note that not all of the Siegel et al. (2007) populations will produce M-giants: Population A in particular is too old and metal poor.¹⁹ All this means, however, is that we would not expect regions of the stellar streams populated by these stars to be traced by M-giants of those metallicities, and we therefore do not calculate K_s magnitudes for these particles.

Zoom In Zoom Out Reset image size

Strictly speaking, the enrichment model presented here is unphysical since it assumes that all of the stellar populations were present in Sgr at all times during its ∼8 Gyr interaction history with the Milky Way, despite the fact that some of these populations did not form until many Gyr later. This is an inescapable consequence of the N-body method adopted in this paper—particles cannot readily be introduced into the simulation at arbitrary times. Ideally, a fully hydrodynamical simulation could build a gas-phase component into the model dwarf, and permit this to form stars at the appropriate times. How to retain an appreciable gas supply in Sgr until only a Gyr ago is a challenging question in its own right, but claims of a detection of gas removed from Sgr as recently as ∼0.4 Gyr ago by its most recent passage through the extended Galactic disk have been made by Putman et al. (2004). Since the N-body particles only represent dynamical tracers, however, we content ourselves with noting that no appreciable quantities of any given stellar population are tidally stripped from the dwarf in this model before they are thought to have formed (see horizontal dashed lines in Figure 14).

Zoom In Zoom Out Reset image size

Tidal stripping of a satellite on an eccentric orbit such as that of Sgr occurs nearly impulsively, with the majority of tidal debris stripped in the brief time interval that the satellite is passing through the pericenter of its orbit. Although the overall trend with successive pericentric passages is to remove stars at successively smaller radii or energies (i.e., "onion-skin" style stripping models, see discussion by Martínez-Delgado et al. 2004), each mass-loss event does not impose a sharp tidal cutoff on the stellar populations of the dwarf. Rather, the rapidly changing tidal radius, the impulsive nature of the stripping, and the dependence of the effectiveness of tidal stripping on the orientation of individual stellar orbits (see, e.g., Read et al. 2006a) means that on any given pericentric passage Galactic tides reach deeply into the satellite and can unbind stars from a variety of populations. This does not completely evacuate the energy/radius range of particles in the satellite, and the remaining stars rapidly relax back to a Plummer model over a dynamical time (∼0.1–0.2 Gyr; see, e.g., discussion by Peñarrubia et al. 2008, 2009). Despite the sharp divisions that we introduced in the distribution of stellar populations within the model satellite, there is significant overlap in the populations expelled from the satellite into the tidal tails over the lifetime of the dwarf. This is illustrated in Figure 14: while particles assigned to Population A are almost entirely stripped from the dwarf during the first two to three orbital periods, stars in Population B are present in many successive epochs of tidal debris and still have a residual bound population (black points at t_unbound = −1) in the dwarf at the present day. Stars in Population D are present in only the most recent tidal debris. The trend of stellar populations along the tidal streams is further affected by the angular phase smearing caused by the differential orbital velocities of stripped stars from each epoch of tidal debris.

By adjusting the relative sizes of each of the stellar populations (i.e., sliding the vertical lines in Figures 13 and 14), it is possible to adjust the relative balance of stars of a given metallicity present in each debris era (i.e., point "color"), and thereby as a function of orbital longitude along the Sgr tidal streams. In Figure 15, we plot the metallicity of stars observed by Chou et al. (2007) and Monaco et al. (2007; their SARG subsample) as a function of orbital longitude, overplotted with the mean metallicity of the simulated Sgr stream using three different assumptions about the energy ranges E* assigned to each of the Populations A–E. The simplest possible assignment in which Populations A/B/C/D/E correspond to uniformly wide energy ranges E* = 1.0–0.8/0.8–0.6/0.6–0.4/0.4–0.2/0.2–0.0, respectively, fails to match the observational data. As illustrated by the red curves in Figure 15, the highest-metallicity debris is overrepresented in the regions of the stream nearest the satellite (i.e., Λ_☉ ≈ 250° and 100° in the leading/trailing streams, respectively), and Population C is overrepresented relative to Populations A and B at larger angular separations (i.e., the model stream is too metal rich at Λ_☉ ≈ 30° and 250° in the leading/trailing streams, respectively). It is possible to suppress the presence of the youngest stellar populations in the Sgr stream by reducing the range in E* assigned to them, and the green curves in Figure 15 show the result of adopting energy ranges E* = 1.0–0.8/0.8–0.6/0.6–0.2/0.2–0.04/0.04–0.0 for Populations A–E, respectively. While this improves the agreement in the regions of the tidal tails nearest Sgr, the mean model metallicity is still too high at larger angular separations.

Zoom In Zoom Out Reset image size

We find that a reasonable fit to both the low- and high-angular separation observational data is achieved with a choice of energy ranges E* = 1.0–0.8/0.8–0.45/0.45–0.38/0.38–0.29/0.29–0.20/0.20–0.04/0.04 − 0.00 for Populations A/B/C1/C2 /C3/D/E, respectively (black lines in Figure 15). Given the large spread of the observational measurements and the large number of assumptions made by our debris-tagging method, a more precise determination of the best-fit choice of energy ranges is not warranted. Rather, we conclude that for some choice of the energy ranges it is possible to broadly fit the observational trends: The mean metallicity of the model Sgr stream decreases from [Fe/H] ∼−0.4 in the core to [Fe/H] ∼−1.0 at large angular separations.

The above choice of energy ranges corresponds to a reasonably strong ∼2 dex metallicity gradient with radius in the initial model Sgr dSph (Figure 16, left-hand panel), consistent with the prediction of Chou et al. (2007), and slightly stronger than observed in the Sculptor dSph by Tolstoy et al. (2004). The strong initial gradient is largely erased by the process of tidal stripping, leaving a relatively weak metallicity gradient ∼0.8 dex remaining by the present day (Figure 16, right-hand panel). The MDF of the bound core of the model Sgr dwarf (Figure 16, right-inset black histogram) is therefore significantly different than both the MDF of the original dSph (Figure 16, left-inset black histogram) and the MDF of the stellar tidal debris that Sgr has contributed to the Galactic halo (Figure 16, right-inset gray histogram). If such strong tidal evolution of the MDF is common for Galactic satellites, it may partially explain the observed chemical abundance mismatch between Galactic halo stars and the present-day dSph population (e.g., Geisler et al. 2005, and references therein). Rather than requiring the early accretion onto the Milky Way of metal poor, high [α/Fe] satellites that are completely disrupted by the present day (e.g., Font et al. 2006), the metal-poor halo stars may simply have originated in the dSphs with which we are familiar, and in which few traces of the originally sizeable metal-poor populations now remain (e.g., Majewski et al. 2002; Muñoz et al. 2006; Chou et al. 2010).

Zoom In Zoom Out Reset image size

One of the challenges in attempting to trace Sgr debris in the 2MASS M-giant view originally provided by Majewski et al. (2003) was that as the metallicity of stars evolves along the stream, so does the zero point of the photometric distance calibration. Since the distance calibration was tied to the relatively metal-rich Sgr core, distances become systematically more inaccurate with increasing separation from the core. Now that we have a model for the evolution of metallicity along the streams, we can perform the comparison of simulations and observational data in a purely observational apparent K_s magnitude space.

In Figure 17 we plot the apparent K_s magnitude, angular position (Λ_☉, B_☉), and radial velocity v_GSR for the simulated Sgr dwarf, overlaid with the Chou et al. (2007) and Monaco et al. (2007) data. We note that the N-body model generally succeeds in reproducing the radial velocity trends indicated by these data, although the Chou et al. (2007) trailing/leading arm sample of stars wrapped into the north/south Galactic hemisphere (crosses/filled squares in Figure 17) may constitute evidence for a more massive Sgr progenitor in order to produce extended wraps of the stream that are broader in radial velocity space. The N-body models also suggest that some of the observational data may originate elsewhere than the Sgr stream: a few of the crosses, three to five filled circles, and two open circles in Figure 17 have radial velocities strongly inconsistent both with the other data points and the N-body models.

Zoom In Zoom Out Reset image size

The agreement between the apparent K_s magnitudes of the N-body model and the Monaco et al. (2007) data appears to be relatively good. Unfortunately, it is difficult to evaluate this agreement for the Chou et al. (2007) data since these authors focused on the brightest available M-giants, which are expected to be biased several or more σ brighter than the mean for the stream (see discussion by Chou et al. 2010). A more instructive comparison is to the M-giant survey of Majewski et al. (2003). In Figure 18, we plot the apparent K_s magnitude for the simulated Sgr dwarf, overlaid with the original M-giant sample from Majewski et al. (2003, i.e., similar to their Figure 8) and the subsample with radial velocity data which we culled from Law et al. (2004) and Majewski et al. (2004) in Figure 1 (open boxes in Figure 18). The K_s magnitudes derived from the model match the larger M-giant sample of the dynamically young sections of the tidal streams (i.e., Λ_☉ ≲ 150° for the trailing arm, and Λ_☉ ≳ 220° for the leading arm) well: over 95% of the radial velocity subsample have K_s magnitudes consistent with the N-body model. The remaining 5% may be Galactic interlopers at smaller distances whose colors and radial velocities happened to be coincident with those of the Sgr stream, although a few located around (Λ_☉, K_s) = (280°, 9.5) may be associated with the T1 wrap of the trailing arm stream instead of the fainter (and more distant) L1 stream since the radial velocity trends of the L1 and T1 wraps are expected to cross in this region of the sky.

Zoom In Zoom Out Reset image size

There are a few regions of stars in Figure 18 that are not explained by the model however. In particular, there is a distribution of M-giants at K_s ∼ 13–14 at all orbital longitudes that are not reproduced by the N-body model. At Λ_☉ ∼ 180°–200° these faint M-giants arguably appear to join smoothly to the trailing arm trend (see also Newberg et al. 2003), suggesting (Pakzad et al. 2004) that the trailing arm may not turn around at this location to approach closer to the Sun again (as the N-body model does) but instead forms a wrap around the entire sky at relatively constant K_s ∼ 13–14. Such a wrap is strongly inconsistent with the predictions of all N-body models to date and, if genuine, may be evidence for a strong evolution in the orbital path of Sgr ∼3 Gyr ago.

7.4.1. Relation to Other Tracers

7.4.2. Uniqueness

Cole et al. (2008) used a maximum likelihood method to characterize the spatial properties of F-turnoff stars in Stripe 82 of the SDSS. They observed tidal debris located at (α, δ, R) = (3137 ± 026, 00, 29.22 ± 0.20 kpc), corresponding to Sgr longitudinal coordinates (Λ_☉, B_☉) = (1052, − 115). As illustrated by Figure 19 (filled circles), the T1 wrap of the N-body model is an extremely good match to these F-turnoff stars in both angular position and distance (matching to within 07 in B_☉, and 1 kpc in distance). Cole et al. (2008) measured the width of the stream perpendicular to the line of sight at this Λ_☉ to have a FWHM =6.74 ± 0.06 kpc. The N-body stream is slightly fatter with FWHM =7.9 kpc, but given the single component mass-follows-light model employed for Sgr the agreement with the apparent width of a specific stellar population is extremely good.

The radial velocity of a sample of RRL in Stripe 82 has recently been measured by Watkins et al. (2009), who detect a component with v_GSR = −130 km s⁻¹ , which is matched by the mean radial velocity of the T1 stream at this position to within 5 km s⁻¹ . Watkins et al. (2009) also measure the metallicity of these RRL, finding that 〈[Fe/H]〉 = −1.41 ± 0.19. This value is appreciably more metal poor than either the M-giant sample or the mean of the N-body model at this longitude ([Fe/H] ∼ −0.7). This difference in measured metallicity emphasizes the strong sensitivity of derived metallicities of a stellar stream to the biases imposed by the type of stellar tracer used (e.g., metal-poor populations do not generally form M-giants while metal-rich populations do not generally form RRL). Thus, until an unbiased assessment of the MDF along the Sgr stream is obtained, models such as that presented in Section 7, which was constrained by a single type of metallicity-biased population tracer, must be considered as merely a demonstration of how MDF variations can originate.

Vivas et al. (2005) published VLT spectroscopy of 16 RRL from the QUEST survey in the range Λ_☉ ∼ 270°–300°. All but 1–2 of these stars, (which have highly discrepant velocities, at v_GSR = −106 and −197 km s⁻¹ ) are in agreement with the angular position, distance, and radial velocity trend of the L1 arm of the N-body model (see Figure 19, open boxes). Starkenburg et al. (2009) have also detected an overdensity of five K-giant stars at similar orbital longitude that are well matched to the distance and radial velocity of the N-body model stream (Figure 19, skeletal triangles), although they appear to be located toward the edge of the stream in B_☉.²⁰

Both samples of stars are well reproduced by the N-body model, although we note that the Vivas et al. (2005) and Starkenburg et al. (2009) samples have mean metallicity 〈[Fe/H]〉 = −1.77 and −1.68, respectively. These metallicities are entirely consistent with each other, but significantly more metal poor than the M-giant sample, which has 〈[Fe/H]〉 ∼ −0.9. As was the case for SDSS Stripe 82, this reinforces the conclusion that stellar populations significantly more metal poor than the M-giant subsample exist at a range of orbital longitudes throughout the Sgr tidal streams. As expected on the basis of dynamical age, the RRL of Vivas et al. (2005) which correspond to tidal debris ∼1.5–3 Gyr old (i.e., magenta points in Figure 19) are slightly more metal poor than the RRL in SDSS Stripe 82 which correspond to ∼0–1.5 Gyr old debris.

We note that an additional sample of metal-poor K-giant stars in a similar longitude range was also pointed out by Kundu et al. (2002; filled triangles in Figure 19). While the L2/T2 wraps of the N-body stream match the radial velocity trend of these stars, the angular positions are discrepant by up to 20° for L2 and 30° for T2. At a typical distance ∼5 kpc, the Kundu et al. (2002) stars are also much closer than expected for the T2 arm (d ∼ 50 kpc), but may be consistent with the L2 arm (d ∼ 10 kpc). If these stars genuinely belong to the Sgr stream, it is most likely that they belong to the L2 wrap, and may suggest that the angular coordinates and/or distance of this wrap in the N-body model require adjustment.

Recent years have witnessed a wealth of substructure discovered in the direction of Virgo, which we follow Martínez-Delgado et al. (2007) in classifying into the Virgo Stellar Stream (VSS) and the Virgo Overdensity (VOD).

The VSS (as traced by QUEST RRL; Duffau et al. 2006) is a relatively well-defined feature that lies at an angular position (Λ_☉, B_☉) = (265°, 14°), i.e., projected on the outer edge of the L1/L2 streams and significantly off the T1/T2 streams in sky position. At a distance of 19 kpc, it is roughly coincident with the L2/T1 stream, but has a radial velocity (v_GSR = 99.8 km s⁻¹ ; Duffau et al. 2006) that is discrepant with these wraps by 220/94 km s⁻¹ , respectively. Curiously, the proper motion of the VSS ([μ_αcosδ, μ_δ] = [ − 4.85 ± 0.85, 0.28 ± 0.85]; Casetti-Dinescu et al. 2009) is reproduced by the T1 stream to within 1σ ([μ_αcosδ, μ_δ] = [ − 4.18 ± 0.41, 0.04 ± 0.25]), but given the angular position and radial velocity mismatch we conclude that the VSS is unlikely to be associated with the Sgr stream.

In contrast, the VOD (Jurić et al. 2008) is a much less well-defined, diffuse clump of stars around (l, b) = (300°, 65°) that covers more than ∼1000 deg² on the sky. The possible association of the VOD with Sgr has been addressed in detail by Martínez-Delgado et al. (2007), who noted that it was roughly coincident with tidal debris from models of the Sgr dwarf (LJM05) disrupting in an axisymmetric, oblate Galactic halo.

In the triaxial halo model derived in this contribution, we find that the model Sgr stream arcs significantly over the solar neighborhood, and tidal debris in the L1 tidal stream at the angular coordinates of the VOD lies much too far away (∼46 kpc for the Sgr stream versus ∼5–17 kpc for the VOD) to be associated (see Figure 19, filled squares). Similarly, both the T1 and T2 wraps lie too far away to produce the VOD, in addition to lying in the opposite angular direction B_☉.

Interestingly, this angular position and distance does correspond closely to the predicted location of the secondary L2 wrap of leading Sgr tidal debris, suggesting that the VOD may be evidence for old tidal debris torn from Sgr ≳3 Gyr ago. If such an identification is correct, we should expect the radial velocity of the VOD to be v_GSR = −131 ± 22 km s⁻¹ . Recent observational work on QUEST RRL in the direction of the VOD by Vivas et al. (2008) derived radial velocity signatures at v_GSR = +215, − 49, − 171 km s⁻¹ . The first (and by far most significant) two of these velocity structures are strongly discrepant with the expected signal by ∼345 and ∼80 km s⁻¹ , respectively. The third structure is possibly consistent with the velocity signature expected of the L2 Sgr stream, differing by only ∼40 km s⁻¹ in this part of the simulated debris where the N-body stream is ill-constrained. Only 7% of the QUEST RRL (i.e., 3/43) are contained in this v_GSR = −171 km s⁻¹ velocity structure, however, and we therefore concur with the conclusion of Vivas et al. (2008) that while Sgr may make some minor contribution to the VOD it is not the dominant origin of the stellar overdensity.

An additional sample of RRL stars from the SEKBO (Moody et al. 2003) survey has recently been studied by Keller et al. (2008), with radial velocities presented by Prior et al. (2009). The Prior et al. (2009) subsample focuses on two groups of RRL which they group into the "VSS region" at α = 12.4–14 hr (i.e., Λ_☉ ≈ 260°–290°) and the "Sgr region" at α = 20–21.5 hr (i.e., Λ_☉ ≈ 10°–40°). These two groups of stars are plotted against our N-body model in Figure 19 (open triangles).

Considering first the angular coordinate plot Λ_☉ versus B_☉, we note that neither of these groups of stars are well centered on the predicted path of the N-body tidal stream, lying at the extreme outer edge of the L1/L2 streams and fully outside the angular path of the T1/T2 streams. This is consistent with Figure 15 of Keller et al. (2008): the plane of the SEKBO survey intersects the Sgr plane at an angle and best overlaps at α ≈ 18 hr (i.e., behind the Galactic center) and α ≈ 5 hr. In the regions traced by the Prior et al. (2009) subsample, we confirm that the RRL tend to lie ∼10°–15° away from the bulk of the Sgr stream. It is perhaps unsurprising then that the distances (typically 16–21 kpc) and radial velocities of the RRL do not obviously match any particular branch of the Sgr stream. While Prior et al. (2009) discussed a possible match between some RRL and trailing Sgr tidal debris in the oblate-halo model of LJM05, no such match obviously exists for the updated model presented here.

We note that the distance (∼20 kpc) and radial velocity range (v_GSR ≈ −200 to +200 km s⁻¹ ) of the Λ_☉ ≈ 10°–40° sample of Prior et al. (2009) is similar to that found for a sample of red clump stars by Majewski et al. (1999; crosses in Figure 19). These red clump stars lie near the angular position of the L1 stream, and LJM05 discussed the possibility that their large spread in radial velocities may indicate the range of velocities occupied by the Sgr stream at this longitude. Given the similarities in distance and radial velocity of the Prior et al. (2009) and Majewski et al. (1999) samples, it is at least likely that they trace the same Galactic substructure. If this substructure is confirmed to be the Sgr stream, it will suggest that the L1 stream at Λ_☉ ∼ 25° has a greater angular width, smaller distance, and larger spread in radial velocities than predicted by the current N-body model.

Keller (2010) has recently described a population of subgiant stars in the range α ≈ 120°–180° from the SDSS that may correspond to the L1 stream of Sgr. While this paper gives insufficient information to include these stars in Figure 19, we note that the Λ_☉ versus distance trend illustrated in Figure 7 of Keller (2010) is reproduced by the revised N-body model to within ∼1 kpc, in contrast to previous N-body models of the Sgr stream (e.g., LJM05). Further characterization of this subgiant population may therefore provide even greater insight into the nature of the Sgr leading tidal stream.