<figure><img src="image_1.png" /></figure>

# Around the Way: Testing ΛCDM with Milky Way Stellar Stream Constraints

Biwei Dai¹,² ID, Brant E. Robertson² ID, and Piero Madau²,³ ID

¹Department of Physics, Peking University, Beijing, 100871, People's Republic of China

²Department of Astronomy & Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

³Institut d'Astrophysique de Paris, Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France

Received 2017 July 21; revised 2018 February 21; accepted 2018 March 31; published 2018 May 9

## Abstract

Recent analyses of the Pal 5 and GD-1 tidal streams suggest that the inner dark matter halo of the Milky Way is close to spherical, in tension with predictions from collisionless N-body simulations of cosmological structure formation. We use the *Eris* simulation to test whether the combination of dissipative physics and hierarchical structure formation can produce Milky Way–like galaxies whose dark matter halos match the tidal stream constraints from the GD-1 and Pal 5 clusters. We use a dynamical model of the simulated *Eris* galaxy to generate many realizations of the GD-1 and Pal 5 tidal streams, marginalize over observational uncertainties in the cluster galactocentric positions and velocities, and compare with the observational constraints. We find that the total density and potential of *Eris* contributed by baryons and dark matter satisfies constraints from the existing Milky Way stellar stream data, as the baryons both round and redistribute the dark matter during the dissipative formation of the galaxy, and provide a centrally concentrated mass distribution that rounds the inner potential. The *Eris* dark matter halo or a spherical Navarro–Frenk–White dark matter work comparably well in modeling the stream data. In contrast, the equivalent dark matter–only *ErisDark* simulation produces a prolate halo that cannot reproduce the observed stream data. The ongoing *Gaia* mission will provide decisive tests of the consistency between ΛCDM and Milky Way streams, and should distinguish between models like *Eris* and more spherical halos.

**Key words:** dark matter – galaxies: fundamental parameters – galaxies: halos – galaxies: kinematics and dynamics – galaxies: spiral – galaxies: structure

## 1. Introduction

The standard ΛCDM model successfully explains a wide range of cosmological observations on scales $k \lesssim 10$ h Mpc$^{-1}$. On smaller scales, comparisons between cosmological N-body simulations and observations have identified several seeming inconsistencies, including the cusp-core (e.g., Flores & Primack 1994; Moore 1994), missing satellite (Klypin et al. 1999; Moore et al. 1999), and Too Big to Fail (Boylan-Kolchin et al. 2011) problems that involve theoretical predictions for the spherically averaged radial density profile or substructure of dark matter halos. These studies have emphasized the need to carefully consider the results of dissipationless N-body simulations when comparing with data, and indeed recent studies with hydrodynamical simulations suggest that baryonic physics may alleviate some of the apparent tensions (e.g., Governato et al. 2012; Pontzen & Governato 2012; Brooks & Zolotov 2014; Madau et al. 2014; Chan et al. 2015; Oñorbe et al. 2015; Wetzel et al. 2016). The shape of the halo represents an important additional prediction of galaxy formation simulations testable through observations. In this paper, we explore whether cosmological simulations of Milky Way–like galaxies can produce realistic dark matter halo shapes when modeling the dissipative effects of baryonic physics.

Collisionless simulations predict triaxial or near-prolate halos, with typical density profile axis ratios [$b/a, c/a$] $\sim 0.4-0.6$ (e.g., Frenk et al. 1988; Dubinski & Carlberg 1991; Warren et al. 1992; Cole & Lacey 1996; Jing & Suto 2002; Vera-Ciro et al. 2011). For external galaxies, observations utilizing gravitational lensing (e.g., Hoekstra et al. 2004; Parker et al. 2007) or X-ray isophotes (e.g., Buote & Canizares 1998; Buote et al. 2002) can test such predictions for halo shapes. These studies find typical halo ellipticities of $\bar{e} \sim 0.35$, in agreement with the results from numerical simulations. Shape

constraints for the Milky Way dark matter halo mostly derive from measurements of the density distribution and kinematics of halo stars or the dynamics of tidal streams such as Sagittarius. Current constraints vary, with studies showing that the Milky Way dark matter halo features an oblate (e.g., Johnston et al. 2005; Loebman et al. 2014), close to spherical (e.g., Ibata et al. 2001; Fellhauer et al. 2006; Smith et al. 2009), prolate (e.g., Helmi 2004; Bowden et al. 2016), or triaxial (Law & Majewski 2010) shape.

Constraints on the shape of the Galactic dark matter halo derived from either stellar halo properties or the Sagittarius stream face complexities from dynamically hot stellar populations, unknown accretion histories, and large angular extents on the sky. In contrast, the comparably dynamically cold tidal streams of the globular clusters Pal 5 (Odenkirchen et al. 2001) and GD-1 (Grillmair & Dionatos 2006) lend themselves more readily to dynamical modeling, and have yielded constraints on the shape of the inner halo of the Milky Way (Koposov et al. 2010; Bowden et al. 2015; Küpper et al. 2015; Pearson et al. 2015; Bovy et al. 2016; Ibata et al. 2016, 2017). Bowden et al. (2015) used a streakline method to model GD-1 in a logarithmic potential model and found a total potential minor-to-major axis ratio of $q_z = 0.90_{-0.10}^{+0.05}$. Küpper et al. (2015) used the same approach to model Pal 5 in a three-component potential model consisting of a bulge, a disk, and a halo, and derived a halo potential minor-to-major axis ratio of $q_{z,h} = 0.95_{-0.12}^{+0.16}$. Using the action-angle method of Bovy (2014), Bovy et al. (2016) modeled the combined data from GD-1 and Pal 5 to determine the axis ratio of the dark matter halo density distribution. Assuming a three-component potential model, they inferred a halo density profile axis ratio of $c/a = 1.05 \pm 0.14$, corresponding to a halo potential flattening close to $q_{z,h} = 1.0$. These results show consistency with
each other but appear in tension with the predictions for halo shapes from collisionless N-body cosmological simulations. Our work revisits the constraints the Pal 5 and GD-1 stellar streams place on the shape of halo in the context of a self-consistent hydrodynamical cosmological model for the formation of Milky Way-like galaxies, the *Eris* simulation (Guedes et al. 2011).

Hydrodynamical simulations have long demonstrated that baryonic dissipation can change the shape of dark matter halos by making them more spherical (e.g., Dubinski 1994; Kazantzidis et al. 2004, 2010; Gustafsson et al. 2006; Debattista et al. 2008; Abadi et al. 2010; Tissera et al. 2010). Dissipation can also redistribute the dark matter mass radially via adiabatic contraction (Blumenthal et al. 1986; Gnedin et al. 2004). The triaxiality of the halo originates in part from box orbits of the dark matter particles. The deepened potential well in the center of a galaxy caused by the condensation of baryons can scatter the box orbits and result in a rounder halo (Dubinski 1994; Debattista et al. 2008). Kazantzidis et al. (2004), Gustafsson et al. (2006), and Abadi et al. (2010) compared the shape of the dark matter halos from collisionless and hydrodynamical simulations, and found that baryonic dissipation causes halos to grow rounder and mildly oblate in their inner regions, with their flattening aligned with any disks formed. Baryonic dissipation increased the halo axis ratios by a few tenths, but did not make them spherical. It therefore remains unclear whether baryonic dissipation can fully explain the discrepancy between halo shapes inferred from stellar streams and the predictions of cosmological N-body simulations. Practical difficulties in forming realistic Milky Way-like galaxies in cosmological hydrodynamical simulations and in comparing the hydrodynamical simulation results with stream data have limited previous studies to simplified models for the dark matter halo structure.

In what follows, we use *Eris* as an example model of a Milky Way-like galaxy and simulate the properties of the Pal 5 and GD-1 streams in this potential. By marginalizing over the uncertainties in the globular cluster positions and velocities, we can determine whether the *Eris* halo shape permits the observed stream orbits. In Section 2, we provide an overview of the *Eris* simulation, provide details on luminous and dark matter properties of the simulated galaxy, and compare these properties with the observed properties of the Milky Way. Section 3 describes our self-consistent field method (SCF), used to efficiently and accurately represent the gravitational potential of *Eris*. In Section 4, we measure and compare the shapes of the dark matter halos from *Eris* and the equivalent dark matter-only simulation *ErisDark* to study the effects of baryonic dissipation. Section 5 describes our method to model and fit the tidal streams of Pal 5 and GD-1 in a variety of potentials, and we present the resulting comparisons with observations in Section 6. We discuss our results in Section 7 and present our conclusions in Section 8.

## 2. The *Eris* Simulation

To generate a cosmologically motivated model for a Milky Way-like potential, we use the results of *Eris*, a cosmological zoom-in simulation of the formation and evolution of a Milky Way-like galaxy (Guedes et al. 2011). The simulation evolved a periodic volume 90 Mpc on a side using the N-body + Smoothed Particle Hydrodynamics code GASOLINE (Wadsley et al. 2004), over the redshift range $z = 90$ to the present in a ΛCDM

<table><caption>**Table 1**
Comparison of the Milky Way and the *Eris* Simulation</caption><thead><tr><th>Property</th><th>Milky Way</th><th><i>Eris</i></th></tr></thead><tbody><tr><td>Dark matter<br>mass <i>M</i><sub>DM</sub> (&lt;20 kpc)</td><td>1.1 ± 0.1 × 10<sup>11</sup> <i>M</i><sub>☉</sub> (a)</td><td>1.1 ± 10<sup>11</sup> <i>M</i><sub>☉</sub></td></tr><tr><td>Total mass <i>M</i><sub>tot</sub> (&lt;60 kpc)</td><td>4.0 ± 0.7 × 10<sup>11</sup> <i>M</i><sub>☉</sub> (b)</td><td>3.34 ± 10<sup>11</sup> <i>M</i><sub>☉</sub></td></tr><tr><td>Disk stellar mass <i>M</i><sub>*,d</sub></td><td>4.6 ± 0.3 × 10<sup>10</sup> <i>M</i><sub>☉</sub> (c)</td><td>3.0 ± 10<sup>10</sup> <i>M</i><sub>☉</sub></td></tr><tr><td>Disk scale length <i>R</i><sub>d</sub></td><td>2.3 ± 0.6 kpc (d)</td><td>2.5 kpc</td></tr><tr><td>Bulge stellar mass <i>M</i><sub>b</sub></td><td>(0.6–1.6) × 10<sup>10</sup> <i>M</i><sub>☉</sub><br>(e), (f)</td><td>0.95 ± 10<sup>10</sup> <i>M</i><sub>☉</sub></td></tr><tr><td>Circular velocity at the<br>Solar circle <i>V</i><sub>c</sub></td><td>218 ± 10 km s<sup>−1</sup> (c)</td><td>212 km s<sup>−1</sup></td></tr><tr><td>Vertical Force at the Solar<br>circle, |<i>z</i>| ≤ 1.1 kpc</td><td>67 ± 6(2πGM<sub>☉</sub> pc<sup>−2</sup>) (c)</td><td>55(2πGM<sub>☉</sub> pc<sup>−2</sup>)</td></tr></tbody></table>

**References.** (a) Bovy et al. (2016), (b) Xue et al. (2008), (c) Bovy & Rix (2013), (d) Hammer et al. (2007), (e) Robin et al. (2012), (f) Portail et al. (2015).

cosmology ($\Omega_M = 0.24$, $\Omega_b = 0.042$, $H_0 = 73$ km s$^{-1}$ Mpc$^{-1}$, $n = 0.96$, $\sigma_8 = 0.76$). The target halo experienced a quiet late merger history and reached a total mass of $M_{vir} = 8 \times 10^{11} M_{☉}$, similar to the Milky Way. The corresponding resampled region contained 13 million high resolution dark matter particles and an equal number of gas particles, for a mass resolution of $m_{DM} = 9.8 \times 10^4 M_{☉}$, $m_{SPH} = 2 \times 10^4 M_{☉}$. The simulation used a gravitational softening length of $\epsilon_G = 120$ pc for all particles, and included Compton cooling, atomic cooling, metallicity-dependent radiative cooling at low temperatures, and heating from a uniform UV background. The star formation algorithm allowed for the stochastic spawning of stellar particles in cold, dense regions, each formed with an initial mass of $m_* = 6 \times 10^3 M_{☉}$. The feedback implementation followed Stinson et al. (2006), with supernova explosions injecting energy and metals in the surrounding gas.

By $z = 0$, *Eris* displays a late-type spiral morphology and resembles a close Milky Way analog. Within the *Eris* virial radius $r_{vir} \approx 235$ kpc lie 7.0 million dark matter, 3.0 million gas, and 8.6 million stellar particles. Table 1 summarizes the basic properties of *Eris* and the corresponding observational constraints of the Milky Way. Most simulated properties show consistency with the observations, assuming a solar galactocentric radius of $R_0 = 8.0$ kpc. The reported measurements of the disk and bulge in the central region follow the orbital circularity-based decomposition by Guedes et al. (2013). The quantity $\epsilon = j_z/j_c$ describes the orbital circularity of the baryonic (gas and stellar) particles, and depends on the angular momentum of each particle along the z-direction $j_z$ and the maximum angular momentum $j_c = \sqrt{GM(<r)r}$. In this decomposition, stellar particles that lie within a radius $r < 2$ kpc and occupy kinematically warm orbits ($\epsilon < 0.5$) get assigned to the bulge, while the remaining particles comprise the stellar disk. For more details about the *Eris* simulation, we refer the reader to Guedes et al. (2011).

For comparison, we also use an equivalent, dark matter-only simulation called *ErisDark* (Pillepich et al. 2014). The detailed properties of the dark matter halos in *Eris* and *ErisDark* are discussed Section 4, and a comparison of their properties enable us quantify the effects of baryonic dissipation. We refer the reader to Pillepich et al. (2014) for other analyses of the dark matter halos formed in the *Eris* and *ErisDark* simulations.
### 3. Methodology

In comparing the Pal 5 and GD-1 stellar stream data to the *Eris* simulation, we must compute the gravitational forces at the positions of stars in the streams. To marginalize over observational uncertainties in the Pal 5 and GD-1 galactocentric positions and velocities, we perform numerical simulations of the stream formation and measure the resulting likelihood of the observed stream position and velocity data, given the *Eris* potential model. Since this approach requires evaluating the simulated halo potential, $\Phi$, many 10s of billions of times, we need an efficient method determining $\Phi$ free from the computation burden of using the simulated particle data. The SCF method (e.g., Clutton-Brock 1973; Hernquist & Ostriker 1992) allows us to construct an accurate model of the *Eris* potential from the particle data using a basis set expansion. We then use this SCF model for the *Eris* potential to simulate the evolution of the stellar streams in a manner analogous to the idealized analytical potentials used by, for example, Pearson et al. (2015) and Bovy et al. (2016). By including contributions from the baryonic components of *Eris*, our approach extends the methodology of Ngan et al. (2015, 2016), who used an SCF model to represent the Via Lactea II (Diemand et al. 2008) halo potential to compute both tidal debris from accreted systems and tidal streams from disrupted globular clusters.

#### 3.1. SCF Method

The SCF method models dynamical systems using a complete set of bi-orthogonal density-potential pair basis functions

$$ \rho(r, \theta, \phi) = \sum_{n,l,m} A_{nlm} \rho_{nlm}(r, \theta, \phi) \quad (1) $$

$$ \Phi(r, \theta, \phi) = \sum_{n,l,m} A_{nlm} \Phi_{nlm}(r, \theta, \phi), \quad (2) $$

where $\rho(r, \theta, \phi)$ represents the density at spherical coordinates $[r, \theta, \phi]$, $\Phi$ is the corresponding gravitational potential, and $A_{nlm}$ are numerical coefficients that weight the contribution of the density ($\rho_{nlm}$) and potential ($\Phi_{nlm}$) eigenfunctions to the expansion.

Representing the density and potential of an arbitrary structure exactly requires an infinite number of terms in the expansion. In practice, one truncates the expansion to remove high frequency terms with high “quantum” numbers $nlm$, resulting in a smoothed potential. To maximize the efficiency of the expansion, it is advantageous to select basis functions such that the first few low-quantum number terms provide a good approximation to the modeled structure. Several analytical basis sets have been proposed (Clutton-Brock 1973; Hernquist & Ostriker 1992; Zhao 1996; Rahmati & Jalali 2009), but the low order radial eigenfunctions of these models do not follow the shape of the Navarro-Frenk-White profile (NFW; Navarro et al. 1996, 1997), with an inner logarithmic slope of $\gamma \sim -1$ transitioning to $\gamma \sim -3$ in the exterior. Following Weinberg (1999), we construct a set of density-potential pair basis functions for the dark matter halo by solving the Sturm-Liouville equation numerically. The basis set functions were tabulated with the radius and stored in lookup tables for interpolation.

To compute the basis functions, we assume their form in spherical coordinates as

$$ \Phi_{nlm}(r, \theta, \phi) = \Phi_{\text{NFW}}(r) u_{nl}(r) Y_{lm}(\theta, \phi) \quad (3) $$

$$ \rho_{nlm}(r, \theta, \phi) = \lambda_{nl} \rho_{\text{NFW}}(r) u_{nl}(r) Y_{lm}(\theta, \phi), \quad (4) $$

where $u_{nl}(r)$ are radial eigenfunctions, $Y_{lm}(\theta, \phi)$ are spherical harmonics, $\lambda_{nl}$ are unknown constants, and $\rho_{\text{NFW}}(r)$ and $\Phi_{\text{NFW}}(r)$ are the NFW density and potential profiles. Each density-potential pair should satisfy the Poisson equation

$$ \nabla^2 \Phi_{nlm} = 4\pi G \rho_{nlm}. \quad (5) $$

This constraint leads to a standard Sturm-Liouville equation,

$$ -\frac{d}{dr} \left[ p(r) \frac{du_{nl}}{dr} \right] + s(r) u_{nl} = \lambda w(r) u_{nl}, \quad (6) $$

where

$$ p(r) = \log^2(1 + r) \quad (7) $$

$$ s(r) = \frac{\log(1 + r)}{(1 + r)^2} + \frac{l(l + 1)}{r^2} \log^2(1 + r) \quad (8) $$

$$ w(r) = \frac{\log(1 + r)}{(1 + r)^2}. \quad (9) $$

With appropriate boundary conditions, these equations were solved to determine $u_{nl}$ and $\lambda_{nl}$ numerically using Matslise 2.0 (Ledoux & Van Daele 2016). The boundary point at $r = 0$ is singular, while the outer boundary is set to $u'_{nl} = 0$ at $r_{\text{max}} = 32r_s$. Since the Milky Way virial concentration ratio is $c_{\text{vir}} \equiv r_{\text{vir}}/r_s \ll 32$, and the region of interest for this calculation (where tidal streams are located) is at galactocentric radii $r \sim r_s$, placing the outer boundary condition at $32r_s$ has little affect on our results. We assume that the potential is symmetric about the $\theta = \pi/2$ plane and azimuthally symmetric. These assumptions allow us to keep only the even $l$ and $m = 0$ terms in the expansion. This approach provides a fast way to compute the potential from the particle distribution. For the *Eris* simulation, this method proves two or three orders of magnitude faster than using a direct summation, depending on which angular symmetries for the potential are assumed. We verified our method by computing expansions for the Hernquist (1990) model and the triaxial NFW model with constant axis ratio

$$ \rho(m) = \frac{\rho_s}{m/r_s(1 + m/r_s)^2}, \quad (10) $$

where

$$ m^2 = x^2 + \frac{y^2}{(b/a)^2} + \frac{z^2}{(c/a)^2}, \quad (11) $$

and $a, b$, and $c$ are the major, intermediate, and minor axis scale lengths. As expected, our SCF method works well when modeling triaxial NFW halos but proves inefficient for highly flattened systems, owing to the use of a spherical NFW profile as the lowest-order term. For modeling the disk component of *Eris*, we include the additional contributions to the potential described below.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Density (left two panels) and potential (right two panels) fields of the dark matter halo in the *Eris* simulation. The first and third panels are calculated from a self-consistent field (SCF) model for the *Eris* potential, while the second and fourth are measured directly from particles. The major axis of the dark halo is aligned with the *x*-axis, and the minor (disk angular momentum) axis is aligned with the *z*-axis. The dashed lines display the best fit elliptical isodensity and isopotential contours calculated from the SCF model.</figcaption></figure>

## 3.2. Eris Potential Model

At least four components are needed to model the complex structure of *Eris*: a dark matter halo, a stellar halo, a bulge, and a baryonic (gaseous and stellar) disk. To perform our decomposition of the *Eris* potential, we adopt the potential minimum computed by the simulation as the center of the galaxy. The total angular momentum vector of all baryonic particles at distances $r < 5$ kpc from the center is calculated, and its orientation is used to define the *z*-axis of the galactic disk. We apply a rotation to diagonalize the inertial tensor constructed from the particle positions in the *x*–*y* plane, completing the coordinate system for our expansions. In constructing the separate bulge, halo, and disk expansions for *Eris*, we aim to provide a quality representation of the total galactic potential, and therefore separate the baryonic particles into distinct structures that approximate the low order terms in each density profile expansion. We first select baryonic particles in an exponential disk distribution and fit their density profile with a combination of three Miyamoto-Nagai disks (Miyamoto & Nagai 1975; Flynn et al. 1996). The remaining baryonic particles are decomposed into a compact bulge and an extended halo resembling an NFW profile, and then both components are modeled using full SCF expansions. This practical method leads to a “disk” versus “bulge” decomposition that differs from the kinematic decompositions reported in Guedes et al. (2013) and Table 1, but provides an accurate representation of the total potential. The dark matter density and potential in *Eris* are separately modeled using a full SCF expansion, as described above in Section 3.1, and the dark matter halo of *ErisDark* is treated in a similar manner.

Figure 1 shows a comparison of the density and potential field of the *Eris* dark matter halo computed from our SCF expansion (the first and the third panels) and the particle distribution (the second and the fourth panels). In the SCF expansion we include terms with quantum numbers $n \le 10$ and $l \le 10$, for a total of 1331 terms. We use direct summation to compute the particle-based potential. The particle density field was smoothed using a spatially adaptive polynomial kernel, with smoothing lengths for each particle computed from the distance to the nearest 32 neighboring particles. The smoothness of the SCF density distribution relative to the particle density field results from truncating the high order terms in the

expansion, but as the isodensity and isopotential contours illustrate, the shape of the halo is well-reproduced. Increasing the number of terms in the expansion does not significantly improve the agreement between the SCF and particle density and potential fields on angular scales relevant to the stream orbits we attempt to model.

One concern about the SCF expansion involves the contributions of substructures, such as dark matter subhalos, to the coefficients of the expansion. Although the angular resolution of the expansion at $l \sim 10$ does not resolve nearby subhalos, the particles in these substructures do contribute power to the expansion coefficients for the $n$ and $l$ terms we retained and therefore could potentially bias the expansion results. For example, if a massive subhalo appears at the major axis of the halo, then the low resolution expansion may give us a more flattened halo. To check the impact of substructures, we compute a test expansion, where we identify and reduce local density maxima in substructures in *Eris* simulation, and then recompute the coefficients on the smoother particle distribution. This procedure removes about 8% of the total particles throughout the whole halo. We randomly select 1 million points at radii $5$ kpc $< r < 30$ kpc, and measure the relative difference of local potential flattening calculated from the two expansions. The local potential flattening is defined as

$$q_{z,h,\text{local}} = \sqrt{\frac{Z F_R}{R F_Z}}, \qquad (12)$$

which equals the axis ratio of the isopotential contour if the latter is a perfect ellipse. Here $Z$ and $R$ correspond to the vertical and radial cylindrical polar coordinates, while $F_R$ and $F_Z$ reflect the projected components of the gravitational force in those directions.

We find that the presence or absence of substructures does not substantially affect the shape of the potential field. The maximum relative difference of the local potential flattening once substructures are removed is only 0.6%. The high order coefficients of the SCF expansion do change after removing the substructures, but these terms have little contribution to the total potential field. The low order coefficients, such as $A_{000}$ and $A_{020}$, are reduced correspondingly because of the loss of mass, but their ratio $A_{020}/A_{000}$, which strongly correlates with the shape of the potential field, remains almost unchanged (the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Minor-to-major axis ratio (left) and intermediate-to-major axis ratio (right) of our potential models as a function of radius. The dark blue color displays the dark matter halo in *Eris* simulation, medium blue shows the total potential in *Eris*, and light blue represents *ErisDark* simulation. The solid lines show the potential flattening, and the dashed lines display the shape of the density field.</figcaption></figure>

relative difference is 0.3%). We therefore conclude that the presence of substructures does not bias the shape of our halo potential model. We also consider contributions of power from the Poisson noise associated with the finite particle number on the expansion coefficients, and for the high resolution *Eris* simulation, the contribution remains small on the radial and angular scales retained in our expansion.

## 4. Density and Potential Shape of *Eris*

Figure 2 shows the minor/major (left panel) and intermediate/major (right panel) axis ratios for the density and potential in the *Eris* (dark blue) and *ErisDark* (light blue) simulations. We fit the isodensity and isopotential contours with ellipsoids, defining the major, intermediate, and minor axes of the ellipsoids relative to the reference frame formed by the galaxy angular momentum and diagonalization of moment of inertia tensor (see Section 3 above). In the x-y and x-z planes, we fit these contours with ellipses every $\Delta r = 1$ kpc over radii $r = 1–100$ kpc with the center $[x_c, y_c]$, the semi-axes $a$ and $b$, and the tilt angle $\tau$ of ellipses allowed to vary. The best fit $x_c, y_c$, and $\tau$ remain close to zero, suggesting that our descriptions for the galaxy centroid and principle axes are accurate.

At radii $r \sim 10–20$ kpc where the tidal streams would reside, the dark matter-only *ErisDark* simulation produces a prolate halo with a minor/major axis ratio $c/a = 0.45$ and an intermediate-major axis ratio $b/a = 0.55$. In contrast, baryonic dissipation greatly increases the sphericity of the dark matter halo formed in the hydrodynamical *Eris* simulation, resulting in an oblate halo that is nearly axisymmetric in the x-y plane with $b/a = 0.98$. Perpendicular to the disk plane, the minor-major axis ratio in the *Eris* simulation increased to $c/a = 0.65$. These results are consistent with those reported by Pillepich et al. (2014), who defined the halo axes following the iterative method described in Kuhlen et al. (2007). In both *Eris* and *ErisDark*, the gravitational potential is much rounder than the density, as expected from the Poisson equation, with a dark matter potential minor-to-major axis ratio of $q_{z,h} \approx 0.85$ in *Eris* and $q_{z,h} \approx 0.67$ in *ErisDark*.

We highlight several features of Figure 2. The axis ratio of the dark matter halo in *Eris* varies with the radius. Within

$r \lesssim 60$ kpc of the center, the halo grows rounder toward the galactic center. Baryonic infall redistributes the dark matter halo mass distribution. After including the contribution of baryons to the total gravitational potential of *Eris*, the potential axis ratios increase at radii $r \gtrsim 10$ kpc. These radii are well beyond the characteristic radii of the baryonic distribution in *Eris*, such that the baryons provide a point source-like contribution to the potential. In the innermost regions of the galaxy, the baryonic disk substantially flattens the total potential and lowers the sphericity of the total potential below that of the dark matter halo potential. The total potential flattening is $q_z \approx 0.88$ in the *Eris* simulation at radii $r \gtrsim 10$ kpc.

## 5. Dynamical Modeling of the Tidal Streams

Armed with an accurate representation of the gravitational potential in the *Eris* simulation, we can model the evolution of tidal streams stripped from globular clusters orbiting in the simulated halo. We then can compare the properties of these model streams with the observed Pal 5 and GD-1 streams, and check for statistical consistency while marginalizing over the uncertainty in the globular cluster positions and velocities in the Milky Way halo.

### 5.1. Streakline Method

To model the tidal streams of Pal 5 and GD-1, we use the streakline method introduced by Küpper et al. (2012). This method was used by Küpper et al. (2015) and Bowden et al. (2015) to constrain the potential of the Milky Way by modeling Pal 5 and GD-1, respectively. As the satellite moves in the potential of the host galaxy, test particles are released near the Lagrange points and their orbits integrated in the combined potential of the galaxy (computed from the SCF method) and satellite (approximated as a point mass). We ignore the mutual gravitational interaction of test particles and the gravitational force on the satellite from the particles. The streakline method proves significantly more computationally efficient than a full N-body simulation, and since the orbits of the stripped stars are assumed to be dynamically independent, the computation is easy to parallelize.
In our model, the stripped stars separate from the progenitor globular clusters near the Lagrange points. We can define the Lagrange radius (King 1962) as

$$r_t = \left( \frac{GM_{\text{sat}}}{\Omega^2 - \partial^2\Phi/\partial r^2} \right)^{1/3}, \quad (13)$$

where $\Omega$ is the angular speed of the satellite in the galaxy potential, $\Phi$ is the gravitational potential of the host galaxy, and $M_{\text{sat}}$ is the mass of the satellite when the stars are released. For Pal 5 we strip stars from the inner Lagrange point to form a leading stream, and from the outer Lagrange point to form a trailing stream. The tracers in the Pal 5 stream models are released every $\Delta t = 8$ Myr, comparable with the $\Delta t = 10$ Myr used by Küpper et al. (2015). For GD-1, following Bowden et al. (2015) we release the particles at a distance $r = 1.2r_t$ from the satellite to prevent a large number of particles from being recaptured by the progenitor cluster. The tracers in the GD-1 stream models are released every $\Delta t = 2$ Myr. By the present day, more than 1000 tracer particles are left in the stream, comparable with the 500–1000 tracer particles in the GD-1 stream simulated by Bowden et al. (2015). Since no progenitor of GD-1 has been observed, we follow Bovy et al. (2016) in assuming the observed stream is leading and only release stars around the inner Lagrange point when modeling GD-1.

The initial velocities of the released stars are determined by

$$v = v_0 + \delta v, \quad (14)$$

where $v_0$ reflects the radial and angular velocity of the satellite (Küpper et al. 2012), and $\delta v$ represents a random component contributed by the velocity dispersion of the satellite. We assume that the internal velocity distribution of the satellite stars follows a Gaussian distribution with zero mean and standard deviation $\sigma_v$. The mass of the satellite declines as the stars are gradually stripped, as represented by a constant mass loss rate $\dot{M}$. Since we only fit the mean orbit of the streams and do not attempt to model the density of stars along the stream, we release the stripped stars at a constant rate. We do not attempt to model variations in the mass loss rate from the globulars that might subtly affect their orbits. Test particles do not correspond to individual stars, and the mass loss rate of the satellite does not scale with the rate at which we release particles. In total, our model has eight free parameters for the progenitor cluster, including its average mass loss rate, present day mass, and six current phase space coordinates. For Pal 5, constraints on the present day progenitor cluster phase space coordinates and the relation between mass and mass loss rate allow us to reduce the number of free parameters to four, including the distance, proper motions, and mass of the cluster.

Starting with the location and velocity of the progenitors at the present, we integrate each progenitor backward by $\Delta t = 4$ Gyr to set the initial conditions. This chosen timescale is long enough for stripped stars to cover the whole Pal 5 and GD-1 streams by the present day, and cautiously exceeds that used by previous analyses (e.g., Bovy et al. 2016). Increasing this timescale would increase the length of the simulated stream, but will have little influence on the stream properties within the extent of the data. We then integrate orbits forward in time using a leapfrog integrator (e.g., Press et al. 1986). We treat any stars that fall within half of the Lagrange radius of

<table><caption>**Table 2**
Parameters of the Spherical Halo Model</caption><thead><tr><th>Component</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td rowspan="2">Halo</td><td>scale length r<sub>s</sub></td><td>17.4 kpc</td></tr><tr><td>mass M<sub>DM</sub>(&lt;20 kpc)</td><td>1.1 × 10<sup>11</sup> M<sub>☉</sub></td></tr><tr><td>Bulge</td><td>total mass M<sub>b</sub></td><td>0.3 × 10<sup>10</sup> M<sub>☉</sub></td></tr><tr><td rowspan="3">Disk</td><td>total mass M<sub>d</sub></td><td>7.0 × 10<sup>10</sup> M<sub>☉</sub></td></tr><tr><td>scale length R<sub>d</sub></td><td>3.0 kpc</td></tr><tr><td>scale height Z<sub>d</sub></td><td>0.28 kpc</td></tr></tbody></table>

their parent globular cluster as “recaptured,” and remove them from the simulation. Once the systems evolve to the present day, we use local regression with a tri-cubic weighting function to calculate the mean track of the stellar orbits along the stream.

Owing to the non-zero velocity dispersion for the stars released from the globular clusters, some particles will move away from the main stream orbit and introduce noise into its determination. Since both the Pal 5 and GD-1 streams are dynamically cold, we reduce the particle noise by reducing the random velocity component in Equation (14). We fix $\delta v = 0$ for Pal 5 (following Küpper et al. 2015) and set $\sigma_v = 0.25$ km s$^{-1}$ for GD-1. We have verified that reducing the random velocity component does not substantially affect the results of our analysis.

## 5.2. Milky Way–like Potential Models

We compare three different models for the gravitational potential of Milky Way–like galaxies in our stellar stream simulations. In the first, we use the SCF expansion of the *Eris* potential described in Section 3.2. We assume the potential to be azimuthally symmetric and mirror symmetric about the disk plane, leaving only 66 terms in our SCF expansion. In our second model, we use an SCF expansion of the *ErisDark* potential. Although the dark matter halo in this simulation is prolate and triaxial, we take an azimuthal average, resulting in an oblate axisymmetric spheroid with a potential flattening of $q_z \approx 0.8$. We make this assumption to match constraints from studies that suggest the axis ratio of the halo in the disk plane is close to unity (Bovy et al. 2014, 2015). As we demonstrate below, the *ErisDark* simulation does not accurately reproduce the stellar stream data and using a prolate model would just further exacerbate this failure. In our third model, we assume a potential similar to the best fit found by Bovy et al. (2016), including a Miyamoto–Nagai disk and a spherical bulge. However, instead of adopting a triaxial NFW dark matter halo with $c/a = 1.05$, we use a spherical NFW halo with the same scale length and normalization. The spherical halo greatly simplifies the calculation and does not have a strong effect on the resulting stream orbits relative to a nearly spherical triaxial ellipsoid. The detailed parameters of this case can be found in Table 2. In the following analysis we treat results from the spherical halo model as a benchmark for evaluating the success of *Eris*, and compare *Eris* and *ErisDark* to study the importance of baryonic physics for altering the potential in which we evolve the stream orbits.

## 5.3. Data

To evaluate the *Eris* potential shape as a template for a Milky Way–like galaxy, we perform an analysis similar to that presented by Bovy et al. (2016). We utilize constraints
<figure><img src="image_4.png" /></figure>

**Figure 3.** Tidal stream data (points with error bars) for the disrupting Pal 5 globular cluster, along with the *Eris* halo (left panels), best fit spherical halo (middle panels), and the *ErisDark* halo (right panels) models. The models attempt to reproduce the path of the tidal stream on the sky by matching the R.A. α and decl. δ of the stars (upper panels) and how their line of sight velocities $V_{\text{los}}$ depend on R.A. (bottom panels). Our Bayesian analysis uses a Monte Carlo approach, varying the line of sight position, proper motions, and mass of Pal 5 as inputs into a streakline model for generating the tidal streams. Tens of thousands of potential stream orbits are considered and compared with the data shown in this figure. The bands show the inner 67% and 95% credibility intervals on stream orbits in the halo potentials. The solid lines present the maximum likelihood stream track in each potential, while the small gray points show a representative subset of the present day test star particles associated with the maximum likelihood parameters. Both *Eris* and the spherical halo allow for orbits that are highly consistent with the Pal 5 data, even though the total potential flattening is quite different for these two models: $q_z = 0.88$ for *Eris* and $q_z = 0.97$ for spherical halo model at the location of Pal 5 cluster. The dark matter-only simulation *ErisDark* does not produce a comparably good model of the Pal 5 data.

from both the Pal 5 and GD-1 stellar streams. For both data sets, to transform between the Celestial coordinates and Galactocentric coordinates, we fix the solar galactocentric radius $R_0 = 8.0$ kpc, the solar vertical position in the Galactic disk $z_0 = 25$ pc, and the solar velocity with respect to the local standard of rest [$V_R, V_T, V_Z$] = [−11.1, 12.24, 7.25] km s⁻¹ (Schönrich et al. 2010).

### 5.3.1. Pal 5

Figure 3 shows the Pal 5 tidal stream data. We use the position of the stream stars on the sky in right ascension (R.A.) α and decl. δ, the six phase-coordinate location of the Pal 5 cluster from Fritz and Kallivayalil (2015), and the line of sight velocity measurement of the Pal 5 cluster from Kuzma et al. (2015) to constrain the model stream orbits. We treat the line of sight velocity of each star as equal to the mean value of the tidal stream at the R.A. of the star, following Bovy et al. (2016). The relation between the mass $M_{\text{cl}}$ and the mass loss rate $\dot{M}$ for Pal 5, $\dot{M}/M_{\text{cl}} = (-0.7 \pm 0.2)$ Gyr⁻¹, is taken from Odenkirchen et al. (2003). We fit both the leading arm and trailing arm of Pal 5 simultaneously.

### 5.3.2. GD-1

Figure 5 shows the GD-1 tidal stream data as points with error bars. The six-dimensional phase space data from Koposov et al. (2010) comprise our constraints, and consist of the positions [$\phi_1, \phi_2$] corresponding to a coordinate system oriented along and perpendicular to the stellar stream, the line of sight distance along the stream, the proper motions along and perpendicular to the stream, and the line of sight velocity of each star along the stream. Similar to Pal 5, the line of sight velocity measurement of each star is treated as the mean value of the stream as a function of position along the sky.

## 5.4. Bayesian Analysis

Given the *Eris* potential, which we treat as fixed according to our SCF model, and a set of parameters θ (some of the six phase space coordinates of the progenitor), we can calculate the mean orbit of the stream following the streakline method described in Section 5.1. We then use a Bayesian analysis to determine the likelihood of the stellar stream data given the model, and locate the maximum likelihood of the progenitor properties that produce streams most consistent with the available data.
<table><caption>**Table 3**
Model Priors for Pal 5</caption><thead><tr><th></th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>Distance</td><td>20 kpc</td><td>24 kpc</td></tr><tr><td>μ<sub>α</sub></td><td>-2.482 mas yr<sup>-1</sup></td><td>-2.110 mas yr<sup>-1</sup></td></tr><tr><td>μ<sub>s</sub></td><td>-2.438 mas yr<sup>-1</sup></td><td>-2.076 mas yr<sup>-1</sup></td></tr><tr><td>Mass</td><td>5 × 10<sup>3</sup> M<sub>☉</sub></td><td>2.5 × 10<sup>4</sup> M<sub>☉</sub></td></tr></tbody></table>

We define our log-likelihood as

$$ \log \mathcal{L}(\theta) = -\sum_{i,j} \frac{[x_{\mathrm{data},ij} - f_{ij}(\theta)]^2}{2\sigma_{ij}^2}, \quad (15) $$

where *i* denotes different data types (e.g., the position or line of sight velocity), *j* runs over the number of data points of the same type, *f*<sub>*ij*</sub>(*θ*) represents the model prediction for the data given parameters *θ*, *x*<sub>data,ij</sub> are the data points, and *σ*<sub>*ij*</sub> are the corresponding uncertainties for each observed star. Since the line of sight velocity data are measurements of individual stars rather than the whole stream, we must account for the velocity dispersion of the stars about the stream. We therefore add a stream velocity dispersion *σ*<sub>sys</sub> to the uncertainty in each star's velocity in quadrature:

$$ \sigma_{ij} = \sqrt{\sigma_{\mathrm{obs},ij}^2 + \sigma_{\mathrm{sys}}^2}. \quad (16) $$

We determine the stream velocity dispersion by fitting a cubic polynomial to the stellar data velocities as a proxy for the stream velocity and then calculating their variance about this fit. We find stream velocities dispersions of $\sigma_{\mathrm{Pal5,sys}} = 2.26 \text{ km s}^{-1}$ and $\sigma_{\mathrm{GD1,sys}} = 15.08 \text{ km s}^{-1}$. To provide the sampling for our analysis, calculate the Bayesian evidence, and compute the posterior distribution of each model, we use *Multinest* (Feroz et al. 2009).

Following Bovy et al. (2016), we fix the celestial coordinates of Pal 5 at $\alpha = 229.018, \delta = -0.124$, and fix the line of sight velocity to $V_{\mathrm{los}} = -58.7 \text{ km s}^{-1}$. We assume a flat prior on the proper motion within the observed uncertainty, as shown in Table 3. The measurements of the distance of the Pal 5 cluster range over 20.9–23.2 kpc (Harris 1996; Vivas & Zinn 2006; Dotter et al. 2011), and we therefore assume a flat prior for the line of sight distance of $D \in [20 \text{ kpc}, 24 \text{ kpc}]$. Mass estimates for Pal 5 vary, ranging from $5 \times 10^3 M_{\odot}$ to $2.2 \times 10^4 M_{\odot}$ (Odenkirchen et al. 2002; Küpper et al. 2015), and we adopt a flat prior of $M_{\mathrm{cl}} \in [5 \times 10^3 M_{\odot}, 2.5 \times 10^4 M_{\odot}]$. For GD-1, since no progenitor has been identified, we assume a flat prior on all six progenitor phase space coordinates, the mass, and the mass loss rate. We first performed a test run with a broad prior, and then reset the prior based on the marginalized one-dimensional posterior distribution of the test run. We provide a complete listing of the model priors for Pal 5 and GD-1 in Tables 3 and 4, respectively.

## 6. Results

In this Section we present the results of the dynamical modeling and Bayesian analysis of the tidal streams Pal 5 and GD-1 in our models for the Milky Way potential.

<table><caption>**Table 4**
Model Priors for GD-1</caption><thead><tr><th></th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>φ<sub>1</sub></td><td>0°</td><td>20°</td></tr><tr><td>φ<sub>2</sub></td><td>-3.5°</td><td>0°</td></tr><tr><td>Distance</td><td>9.5 kpc</td><td>15.5 kpc</td></tr><tr><td>μ<sub>φ1</sub></td><td>-9 mas yr<sup>-1</sup></td><td>-4 mas yr<sup>-1</sup></td></tr><tr><td>μ<sub>φ2</sub></td><td>-2.2 mas yr<sup>-1</sup></td><td>-1.2 mas yr<sup>-1</sup></td></tr><tr><td>V<sub>los</sub></td><td>-380 km s<sup>-1</sup></td><td>-250 km s<sup>-1</sup></td></tr><tr><td>Mass</td><td>10<sup>4</sup> M<sub>☉</sub></td><td>10<sup>5</sup> M<sub>☉</sub></td></tr><tr><td>Mass loss rate</td><td>0 M<sub>☉</sub> Myr<sup>-1</sup></td><td>20 M<sub>☉</sub> Myr<sup>-1</sup></td></tr></tbody></table>

<table><caption>**Table 5**
Log-evidence for Milky Way-like Potential Models</caption><thead><tr><th></th><th>Pal 5</th><th>GD-1</th></tr></thead><tbody><tr><td>Spherical Halo</td><td>-43.30 ± 0.09</td><td>-35.08 ± 0.11</td></tr><tr><td>Eris</td><td>-42.20 ± 0.09</td><td>-36.10 ± 0.11</td></tr><tr><td>ErisDark</td><td>-81.39 ± 0.11</td><td>-49.90 ± 0.14</td></tr></tbody></table>

<table><caption>**Table 6**
Best χ² for Milky Way-like Potential Models</caption><thead><tr><th></th><th>Pal 5</th><th>GD-1</th></tr></thead><tbody><tr><td>Spherical Halo</td><td>31.57</td><td>19.01</td></tr><tr><td>Eris</td><td>31.47</td><td>19.11</td></tr><tr><td>ErisDark</td><td>65.77</td><td>25.17</td></tr></tbody></table>

### 6.1. Pal 5

Figure 3 shows the 68% and 95% credibility intervals for the mean tidal stream orbits computed in the *Eris* (left panels), spherical halo (middle panels), and *ErisDark* (right panels) models for the decl. $\delta$ (upper panels) and line of sight velocity $V_{\mathrm{los}}$ (bottom panels) as a function of R.A. $\alpha$. To compute the credibility intervals at each position along the stream, likelihood samples from the Bayesian analysis were used to compute the posterior predictive distributions for decl. $\delta$ or line of sight velocity $V_{\mathrm{los}}$ as a function of R.A. $\alpha$. In calculating the posterior predictive distribution, we assigned uncertainties to regions between the data points by interpolation. As Figure 3 demonstrates, both the *Eris* and spherical halo models allow for stellar streams that appear highly consistent with the available data. We find this compatibility remarkable given that the spherical model adopts the best fit model of Bovy et al. (2016), with a halo mass and shape selected from a full Bayesian analysis to provide maximal consistency with these data, whereas the *Eris* potential was in no way selected to reproduce the Pal 5 stream.

To quantify the relative agreement, we present the log-evidence of both the *Eris* and spherical halo models for the Pal 5 stream in Table 5 and the $\chi^2$ of the best fit parameters in Table 6. The Bayesian evidence is defined as the average of the likelihood over the prior, and can be written as

$$ \mathcal{Z} = \int \mathcal{L}(\theta) \pi(\theta) d\theta, \quad (17) $$

where $\theta$ are the parameters in the model, $\mathcal{L}(\theta)$ is the likelihood, and $\pi(\theta)$ is the prior. *Eris* provides the best model, but the spherical halo model allows for streams that show comparable consistency with the data. Given that the spherical halo model
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Marginalized 1D posterior probability distributions of the three Pal 5 cluster parameters. The priors on these parameters are set using observational constraints (see Section 5.4 and Table 3). The light blue color displays the normalized posterior distribution for the *Eris* potential, while the dark blue color shows the spherical halo model results.</figcaption></figure>

parameters remained fixed during our analysis, this comparison is formally correct, but we note that the additional flexibility of allowing the halo mass and radius to vary could produce a somewhat better log-evidence than shown in Table 5. Nonetheless, Figure 3 and Table 5 suggest that the *Eris* and spherical halo models work comparably well. In Table 5, we also present the log-evidence for models generated from the dark matter–only *ErisDark* simulation. *ErisDark* performs poorly, which is unsurprising given the strongly aspherical potential without dissipative effects or a large central baryonic mass concentration. The end of the leading arm and the trailing arm of the best fit *ErisDark* stream both deviate from the observed data. Stars that lie at the end of the trailing are typically the first stripped and travel the farthest away from the cluster, and therefore are quite sensitive to the shape of the host galaxy potential.

Figure 4 presents the one-dimensional marginalized posterior distribution of the Pal 5 distance and proper motions in the *Eris* and the spherical halo models. As expected, these parameters have quite different posteriors in the *Eris* and spherical halo models, reflecting differences in their potential shapes. The overlap between the two posterior distributions is comparable to the width of the spherical halo posterior with $\sigma_{\mu} \approx 100~\mu\text{as yr}^{-1}$, suggesting that future data more precise than this level would be required to distinguish between the two models. We note that owing to the strong correlation between distance and proper motions, the posterior distribution for the distance to Pal 5 cluster in the spherical halo model truncates at $D \approx 21$ kpc, corresponding to a proper motion near the lower limit of the proper motion prior. However, we have checked that increasing the proper motion prior range does not substantially change these results.

## 6.2. GD-1

Figure 5 shows the 68% and 95% credibility intervals for tidal stream orbits computed in the *Eris* (left panels), spherical halo (middle panels), and *ErisDark* (right panels) models. The data for the GD-1 stream are presented in the rotated spherical coordinates ($\phi_1, \phi_2$), oriented parallel and perpendicular to the direction of the stream orbit. Shown are the positions $\phi_2$, line of sight distance, proper motions, and line of sight velocities as a function of the $\phi_1$ coordinate for each model. Only a representative subsample of 500 tracers are plotted, while

more than 1000 tracers are used to calculate the stream centerline. As with Pal 5, we show the marginalized 68% and 95% credibility intervals for each orbital parameter as bands, and the maximum likelihood orbit as a solid line. Both the *Eris* and spherical halo potential models provide comparable fits, given their complexity. In contrast, the best fit stream in the *ErisDark* potential encounters difficulty at large $\phi_1$ and produces broader streams than the other potential models. Tables 5 and 6 provide the log-evidence and best $\chi^2$ for each model, showing that the spherical halo is slightly favored to *Eris*, but that the difference in log-evidence is less than for Pal 5. The *ErisDark* model performs substantially poorer.

## 7. Discussion

Our analysis demonstrates that the gravitational potential of the simulated *Eris* galaxy permits the formation of tidal streams that closely resemble Pal 5 and GD-1. We emphasize that we did not select *Eris* among a range of models, nor did we tune its physical prescriptions for the purpose of producing a simulated halo that would allow for realistic stellar streams. *Eris* simply provides a close Milky Way analog that serves as a model gravitational potential for comparing with observed Milky Way stellar streams. We did not attempt to find a statistically optimal model potential in the maximum likelihood sense, as the *Eris* potential remained fixed throughout our tests. As we show below, with the greater flexibility of allowing the halo shape or disk mass to vary, we can find a somewhat closer match to either the Pal 5 or GD-1 stellar stream data than *Eris* provides. However, *Eris*, a simulated galaxy formed fully self-consistently in a complete cosmological context, works very well as a Milky Way–like potential model, as probed by the Pal 5 and GD-1 streams. Importantly, this result suggests that structure formation in the $\Lambda$CDM cosmology can produce galaxies whose halos closely resemble that of the Milky Way. The $\Lambda$CDM structure formation process can accomplish this success without resorting to additional physics beyond those associated with standard CDM and the dissipational processes associated with baryons.

### 7.1. Effects of Baryonic Dissipation

The comparison of the dark matter halos in the *Eris* and *ErisDark* simulations in Section 4 demonstrated that baryonic dissipation greatly increased the sphericity of the halo.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Comparison of the average radial density profile ($\rho(r)$; left) and interior mass distribution ($M(<r)$; right) of dark matter halo models used in computing stream orbits. In the right panel, the interior mass distribution of each model is normalized by the NFW interior mass profile. Relative to a spherical NFW halo (lightest blue), baryonic dissipation increases the central concentration of the *Eris* halo mass distribution (darkest blue). As expected, *ErisDark* (dark blue) and a flattened NFW halo (light blue) with $q_{z,h} = 0.85$ (chosen to have oblateness similar to *Eris*) are slightly more centrally concentrated relative to a spherical NFW.</figcaption></figure>

Figure 6 shows the average radial density profile $\rho(r)$ and interior mass distribution $M(<r)$ for *Eris*, *ErisDark*, and a spherical NFW. The prolate halo in the dissipationless *ErisDark* simulation is slightly more centrally concentrated than a spherical NFW halo. The much poorer performance of *ErisDark* in reproducing the streams relative to the spherical NFW halo therefore results primarily from the flattened *ErisDark* halo shape, rather than large differences in their radial mass distribution. In contrast, baryonic dissipation redistributes the mass of *Eris* inward compared with a spherical NFW halo and rounds the halo relative to *ErisDark*.

Since *Eris* performs comparably to a spherical NFW halo but has a flatness intermediate between *ErisDark* and spherical, we performed additional tests to explore the relative roles of radial mass re-distribution or rounding of the halo shape in improving the *Eris* performance. We revisited our previous analysis using a spheroidal NFW halo with a flattening $q_{z,h} = 0.85$. This model mimics the *Eris* halo shape, but uses the more extended NFW halo mass distribution. In terms of the mass $M_s$ within the scale radius $r_s$, the potential of a spheroidal NFW halo can be written as

$$ \Phi^h(m) = -\frac{GM_s}{m} \log(1 + m/r_s), \qquad (18) $$

where

$$ m^2 = x^2 + y^2 + \frac{z^2}{q_{z,h}^2}. \qquad (19) $$

If $q_{z,h} = 1$, then Equation (18) reduces to the spherical NFW potential. The mass distribution of this flattened NFW model is also plotted in Figure 6, showing as expected that a spheroidal NFW model is slightly more centrally concentrated than a spherical NFW model.

Using this flattened NFW potential model with $q_{z,h} = 0.85$, we perform the same Bayesian analysis described in Section 5.4. For the Pal 5 stream, we find a log-evidence of $\log Z = -41.55 \pm 0.09$, slightly better than for *Eris* or a spherical NFW halo. For the GD-1 stream, we find a log-evidence of $\log Z = -38.17 \pm 0.12$, worse than for *Eris* or a spherical NFW halo. The performance of the spheroidal NFW

halo when combining both the Pal 5 and GD-1 constraints is worse than for either *Eris* or a spherical NFW halo. From these results we infer that the mass re-distribution by baryonic infall is somewhat more influential than the rounding of the halo shape in improving the performance of *Eris*, but that the relative importance of these two effects may vary with the radius.

We also note that the degree of dissipation in *Eris* should be model-dependent and in principle the effect of baryonic dissipation could be stronger with, for example, a different feedback prescription or in a galaxy with a different formation time than *Eris*. Kazantzidis et al. (2010) used simulations to show that the amount of sphericalization not only depends on the mass of the forming disk, but also on its alignment with halo. The effect of dissipation in rounding the halo is strongest if the symmetry axis of the disk is aligned with the major axis of the dark halo. Since the axis of the disk in *Eris* is aligned with the minor axis of the halo, and the *Eris* disk mass is smaller than the Milky Way, it is plausible that the effect of baryonic dissipation on the halo sphericity could be stronger in the real Milky Way.

## 7.2. Comparison to Previous Studies

Several previous analyses studied the development of Milky Way stellar streams, and used observed streams as constraints on the Milky Way halo shape. Our formalism enables us to compare straightforwardly with previous analyses and, in some cases, attempt a reproduction of their results. In all cases, we find broad consistency where we can compare. This agreement helps validate both the findings of previous studies and the Bayesian methodology presented in this work.

### 7.2.1. Halo Potential versus Density Shape

In what follows, we compare with studies that constrain either the halo potential flattening $q_{z,h}$ or the halo density flattening, expressed in terms of the minor-to-major axis ratio $c/a$. To make these comparisons clearer, we link the shape of the density field to potential flattening by calculating the potential field of a spheroidal NFW density profile. We model the potentials of flattened NFW profiles with different axis ratios with our basis
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Potential flattening vs. axis ratio of the density field for a flattened NFW profile. The dark blue color shows our measured relation, in which $q_{z,h}$ is measured by fitting ellipse to the potential field at the $r_s$ radius, and the potential field is calculated using our basis expansion. The light blue color displays our best fit line in the form of Equation (20). This relation is only satisfied near the inner halo ($r \sim r_s$) of a NFW profile with a constant density axis ratio.</figcaption></figure>

expansion, and then fit ellipses to the isopotential contours to measure the potential flattening. Figure 7 shows this simple correspondence between $q_{z,h}$ and $c/a$. We find that a quadratic function can fit well this relation as

$$q_{z,h} = 1 + \eta(c/a - 1) - \xi(c/a - 1)^2, \quad (20)$$

with the numerical constants $\eta \approx 0.417$ and $\xi \approx 0.056$. Note that we measured the potential shape at the NFW profile scale radius $r = r_s$, and since the potential flattening is radius-dependent, this relation remains only satisfied in the inner halo with $r \sim r_s$. As expected, a spherical density profile results in a spherical potential. Note that Figure 7 and Equation (20) are consistent with the connection between $q_{z,h}$ and $c/a$ in the *Eris* halo, as shown in Figure 2. Equation (20) also provides a rule of thumb for comparing the relative statistical uncertainties $\sigma_q$ or $\sigma_{c/a}$ on $q_{z,h}$ or $c/a$, and via the error propagation formula, we have

$$\sigma_{c/a} \approx \frac{\sigma_q}{\sqrt{\eta^2 - 4\xi(q_{z,h} - 1)}}. \quad (21)$$

Near $q_{z,h} \approx c/a \approx 1$, we have roughly that $\sigma_{c/a}/\sigma_q \approx \eta^{-1} \approx 2.4$ for NFW halos.

### 7.2.2. Comparison to Previous Streakline Models

The *Eris* simulation produced a Milky Way–like galaxy with an oblate dark matter halo and a total potential flattening $q_z = 0.88$, and our results suggest that the *Eris* potential can produce tidal streams remarkably consistent with the Pal 5 and GD-1 data. This result appears compatible with results from previous analyses using streakline models to simulate the stellar stream orbits. Bowden et al. (2015) used a similar

streakline method to model GD-1 in a logarithmic potential and found a total potential flattening for the Milky Way of $q_z = 0.90^{+0.05}_{-0.09}$, inferring that the Milky Way halo is flattened at the orbital radius of GD-1. Küpper et al. (2015) used a streakline method to model Pal 5 in a three-component potential consisting of a variable NFW halo and a fixed disk and bulge. Küpper et al. (2015) fit the main stream orbit positions, line of sight velocities, and overdensities of the Pal 5 stream to the data. They found that the shape of the Milky Way halo potential is constrained to be $q_{z,h} = 0.95^{+0.16}_{-0.12}$, and conclude the Milky Way halo is spherical. We note that given their quoted uncertainties, both the Bowden et al. (2015) and Küpper et al. (2015) analyses are consistent with the *Eris* halo potential flattening ($q_{z,h} = 0.85$) to within $1\sigma$ and permit quite oblate halos.

### 7.2.3. Comparison to Previous Action-angle Models

Using the Bovy (2014) action-angle method for evolving the stellar stream orbits, Bovy et al. (2016) combined the data from Pal 5, GD-1, and other Milky Way observations to constrain the shape of the Milky Way halo, assuming an NFW density profile. Their Milky Way potential model consisted of a flattened dark matter halo, a Miyamoto–Nagai disk, and a spherical bulge. They allowed for a variable NFW halo profile shape, leaving the halo density axis ratio $c/a$ as a free parameter of their model.

When modeling Pal 5, Bovy et al. (2016) allowed the distance, proper motion, and velocity dispersion of the cluster to vary. We have adopted their priors on the Pal 5 progenitor distance and proper motions (see Table 3 for numerical values). Given the data and priors, Bovy et al. (2016) constrained the halo axis ratio to $c/a = 0.89^{+0.20}_{-0.23}$. They further measured the force at the location of the stream progenitor, and from this inferred that the total potential flattening was $q_z = 0.94 \pm 0.05$. Our *Eris* halo axis ratio and total potential flattening lie at the $1-\sigma$ limit of their Pal 5 constraint.

When modeling GD-1, Bovy et al. (2016) fixed the progenitor cluster internal velocity dispersion to $\sigma_v = 0.4$ km s$^{-1}$, adopted $\phi_1 = 0$ assuming that the stream represents a leading arm, and varied the other five phase space coordinates of the GD-1 progenitor as free parameters. Following a similar analysis to Pal 5, Bovy et al. (2016) constrained the halo density axis ratio to $c/a = 1.27^{+0.27}_{-0.22}$ and the total potential flattening to $q_z = 0.95 \pm 0.04$ for an NFW profile halo.

Bovy et al. (2016) combined the Pal 5 and GD-1 force measurements with other Milky Way potential constraints, and found a combined constraint on the axis ratio of the dark matter halo of the Milky Way of $c/a = 1.05 \pm 0.14$, assuming the halo density profile was NFW. The dark matter halo density axis ratio for *Eris* ($c/a \approx 0.65$) lies $\sim 3\sigma$ discrepant with this result, even as *Eris* and a spherical NFW dark matter halo perform comparably well in reproducing the stellar stream data. We now turn to an investigation of the origin of this apparent tension.

First, we review differences in the baryonic components of each model. Comparing with the best fit model of Bovy et al. (2016), *Eris* has a larger bulge and a smaller disk (see Tables 1 and 2). Therefore the contribution of the baryons to the total potential in the *Eris* model is more spherical and may partially compensate for the non-sphericity of the dark matter halo. To test this hypothesis, we combine the dark matter halo from *Eris* and the baryonic component from the spherical halo model to
form a new Milky Way potential, and calculate the Bayesian evidence of this model to test for consistency with the tidal stream data following the procedure described in Section 5. We find that this “massive disk” model, with a disk mass increased to $M_d = 7.0 \times 10^{10} M_\odot$, performs comparably well as the *Eris* or spherical NFW halo models, with a log-evidence $\log Z = -42.05 \pm 0.09$ for Pal 5 ($\chi^2 = 28.96$) and $\log \bar{Z} = -35.19 \pm 0.12$ for GD-1 ($\chi^2 = 18.88$), comparable to *Eris*. We therefore conclude that the relative mass of the disk and bulge in *Eris* does not strongly affect its ability to reproduce the stellar stream data.

The above results arise due to several physical effects. The differences in the total potential resulting from adopting the *Eris* disk or the “massive disk” adopted from Bovy et al. (2016) are not large at the location of the Pal 5 stellar stream. The potential at the location of Pal 5 is dominated by the contribution from the dark halo, and the influence of the disk there is weak. As we discussed in Section 4, the characteristic length of the disk ($\approx 2.5$ kpc) is much smaller than the radial location of the streams. At $r \approx 15$ kpc, the compact baryonic component acts like a point source and makes the total potential more spherical than the dark matter potential. Therefore increasing the disk mass will induce a more spherical total potential, and using the *Eris* disk in our analysis is a conservative assumption.

We emphasize again that the disk of *Eris* is smaller than the Milky Way disk, and therefore *Eris* does not match all constraints on the Milky Way disk mass in the solar neighborhood. However, we do not argue that *Eris* is a perfect Milky Way model, but instead are trying to understand the robustness of current stellar streams on the shape of the Milky Way halo. Clearly the (dis-)agreement between *Eris* and the Milky Way at the solar circle does not necessarily equate to their (dis-)agreement at 20 kpc. Replacing the *Eris* disk with a more massive one, as we did above, will just serve to make *Eris* a “better” (i.e., more spherical potential) model for a Milky Way–like galaxy that can satisfy the stream constraints at both 8 and 20 kpc.

### 7.2.4. Variable Halo Shape

We expect that the action-angle method (Bovy 2014) and our streakline method produce model streams of comparable physical accuracy, and therefore our method should reproduce the results of Bovy et al. (2016), given similar assumptions. If not, then discrepancy between the *Eris* and spherical NFW halo results could owe to the differences in modeling the streams. We therefore revisited our analysis, allowing for an additional free parameter $q_{z,h}$ describing the potential shape of a spheroidal NFW halo model (see Equation (18)). This analysis generalizes the results discussed in Section 7.1 for a spheroidal NFW halo potential to a variable flattening $q_{z,h}$. We fix the baryonic component, the scale mass, and the scale length of the halo to the best fit value in Bovy et al. (2016), only allowing the potential flattening to vary within $q_{z,h} = [0.7, 1.4]$. Other details of the stream modeling follow the method described in Section 5. Since both Pal 5 and GD-1 lie close to the scale length of the Milky Way halo ($r_s = 18.0 \pm 7.5$ kpc; Bovy et al. 2016), we use Equation (20) to transform the axis ratio constraints in Bovy et al. (2016) to potential shapes for a direct comparison.

We find that Pal 5 constrains the NFW potential shape to $q_{z,h} = 0.87 \pm 0.06$, more flattened but consistent with the

value $q_{z,h} = 0.95_{-0.10}^{+0.09}$ (transformed from $c/a = 0.89_{-0.23}^{+0.20}$) in Bovy et al. (2016). *Eris* lies at the $1\sigma$ limit of their result, but close to the center of our constraint. The difference in these two constraints may owe to our fixing of the baryonic galaxy parameters, while Bovy et al. (2016) allowed for a wide range of both disk and halo model parameters. Overall, we conclude our model can produce results consistent with Bovy et al. (2016) under a similar set of assumptions.

For GD-1 we find that the data prefers a prolate NFW halo with a potential flattening $q_{z,h} = 1.21 \pm 0.11$, in agreement with the value $q_{z,h} = 1.11_{-0.09}^{+0.10}$ (transformed from $c/a = 1.27_{-0.22}^{+0.27}$) found by Bovy et al. (2016). That the GD-1 data alone favor a prolate NFW halo (without additionally considering Pal 5) does not necessarily mean that an oblate halo will fail to reproduce the GD-1 stream data. Even though *Eris* has a smaller evidence as a model for GD-1 than a spherical NFW halo, we can also find a given set of initial conditions for the GD-1 progenitor that would allow the *Eris* potential to reproduce the observed streams with high fidelity (see Figure 5 and Table 6).

Finally, we combine the Pal 5 and GD-1 data and fit them simultaneously using an oblate NFW halo and 13 free parameters (4 from Pal 5, 8 from GD-1, and $q_{z,h}$). We find that the NFW halo potential shape is constrained to be $q_{z,h} = 0.97 \pm 0.07$, consistent with the results from Bovy et al. (2016). This result is also consistent with a simple variance-weighted average of the constraints from Pal 5 and GD-1 separately. We have shown that the oblate *Eris* halo is consistent with the data, but we agree with Bovy et al. (2016) that for an NFW profile, the data favor a spherical halo. Since we have used the same stream modeling method and the same parameter setup, the only difference between these two analyses is the form of the potential (i.e., *Eris* differs from an analytical oblate NFW halo). *Eris* is a dark halo formed self-consistently in a cosmological simulation that models the dissipative physics relevant for the formation of a realistic Milky Way–like galaxy. Given the comparable quality of the fits and the current constraining power of the data, we suggest that both the oblate *Eris* halo and a spherical NFW halo provide reasonable models for dark matter halos in galaxies like the Milky Way.

### 7.2.5. Comparison to Other Tidal Stream Analyses

While we use the *Eris* hydrodynamical simulation to provide the background potential for our stellar stream modeling, prior work has used collisionless dark matter–only simulations to model the development of Milky Way stellar streams. Ngan et al. (2015, 2016) simulated tidal streams in the potential of the Via Lactea II (VL2) N-body simulation and a spherical analytical potential, with and without subhalos, to study the stream morphologies. They used a similar SCF approach to ours for representing the simulated dark matter halo, except they used the Hernquist (1990) radial basis functions instead of the NFW model we adopt. Ngan et al. (2015, 2016) found that orbiting subhalos and the time-dependence of the potential shape could lead to both fluctuations in the density and eventually the dispersal of the tidal streams, consistent with previous studies (e.g., Ibata et al. 2002). They also found that massive subhalos could increase the stream density by shepherding stream stars.

More recently, Sandford et al. (2017) used the VL2 potential and other, more idealized smooth potentials to further examine the interplay between stream morphologies, the global halo
shape, and the presence of substructures. Sandford et al. (2017) found that the shape of the halo potential and substructure could affect stream morphological properties such as their thickness or the regularity of their stellar distribution.

Our work largely complements these previous efforts, in that we only use the positions and velocities of the streams to evaluate *Eris* as a model of a Milky Way–like halo and do not attempt to use the morphological properties of the streams to evaluate the model potential. We do find that the aspherical halo in the dark matter–only *ErisDark* simulation widens the model streams. The end of the Pal 5 trailing arm becomes relatively broad, even though we set the internal velocity dispersion of the progenitor cluster to $\sigma_v = 0$, and the test particles take twice as long cover the whole Pal 5 stream in *ErisDark* than in *Eris*. With a non-zero progenitor internal velocity dispersion, the stream-fanning in *ErisDark* would be exacerbated (Pearson et al. 2015). The same phenomenon also occurs in the GD-1 model.

### 7.3. Gaia and Future Observations

Given the current quality of data, we have demonstrated that both *Eris* and a spherical halo model prove consistent with the current Pal 5 and GD-1 stream data, even though the models differ structurally. Relative to our spherical halo model, *Eris* has a smaller disk, less mass within 20 kpc, and a substantially flattened halo ($q_{z,h} = 0.85$), but both models perform comparably well. However, with more precise data, one or both models might be falsifiable. Based on our fitting to the Pal 5 stellar stream, we constrain the Pal 5 globular cluster parameters to be $[D, \mu_\alpha, \mu_\delta, M] = [21.0 \pm 0.8 \text{ kpc}, -2.337 \pm 0.096 \text{ mas yr}^{-1}, -2.302 \pm 0.100 \text{ mas yr}^{-1}, 1.3 \pm 0.6 \times 10^4 M_\odot]$ in *Eris* and $[21.7 \pm 0.5 \text{ kpc}, -2.393 \pm 0.058 \text{ mas yr}^{-1}, -2.349 \pm 0.061 \text{ mas yr}^{-1}, 1.4 \pm 0.4 \times 10^4 M_\odot]$ in the spherical halo model. We find that proper motion is the most sensitive parameter for distinguishing between different potentials, as further indicated by the one-dimensional marginalized posterior distributions shown in Figure 4. Based on the posterior distribution from the *Eris* model, we estimate that improving the proper motion constraints to $\sigma_\mu \lesssim 50 \mu\text{as yr}^{-1}$ for the streams would distinguish between *Eris* and a more spherical halo, well beyond the current $\sigma_\mu \approx 150 \mu\text{as yr}^{-1}$ uncertainties. According to de Bruijne (2012), *Gaia* should reach final mission proper motion constraints of $\sigma_\mu \lesssim 100 \mu\text{as yr}^{-1}$ for individual stars with the apparent brightness of stars in the Pal 5 and GD-1 streams ($m \approx 17-20$), and will provide much tighter constraints by combining measurements from multiple stars in each stream. We therefore expect the final mission *Gaia* data will be able to distinguish between models like *Eris* and a more spherical halo model for the Milky Way.

### 7.4. Time Evolution of the Eris Halo

Dark matter halos evolve over time owing to accretion. If the shape of the host halo potential changes dramatically during the period of the formation of a tidal stream, using a static halo potential to model the tidal stream could be inappropriate. We modeled the potential of *Eris* halo at redshift $z=0$ using our basis expansion, and we assumed that the halo did not evolve over the formation time of the stellar streams. To test the accuracy of this assumption, we measure the shape and inner mass at $r \lesssim 20 \text{ kpc}$ of the *Eris* halo as a function of time, as shown in Figure 8. The shape of the *Eris* halo has evolved very

slowly over the last 8 billion years. Over the same period, the inner halo increased its mass by only about 10% and remained quite stable over the last 4 billion years. In our model, the Pal 5 and GD-1 stellar streams developed over the last ~4 billion years, and these measurements validate our use of a static halo when using Pal 5 and GD-1 streams to analyze the constraints on the Milky Way potential.

### 7.5. Caveats

The results of our analysis are subject to important caveats. The comparison of simulated galaxies like *Eris* with real systems like the Milky Way often involve an imperfect statistical rationale, as the expense of high resolution cosmological simulations typically limit the range of model galaxies available to compare with observations. Our analysis deals with only a single realistic Milky Way analog, and it remains unclear whether the results we present for *Eris* generalize to all simulated galaxies that resemble the Milky Way. We plan to revisit our results using the *Venus* simulation (see Sokolowska et al. 2017), a similar simulation to *Eris* but where the resimulated galaxy is chosen to have a more violent merger history. Even given these caveats, our results do indicate both the halo mass distribution and its flattening matter in shaping stellar streams, and appropriate caution should be employed when assuming a halo density profile while attempting to infer the flattening of the halo potential.

## 8. Conclusions

We have used the *Eris* hydrodynamical simulation of a Milky Way–like galaxy forming in a full cosmological context (Guedes et al. 2011) as a model for reproducing the observed tidal streams Pal 5 and GD-1. We find that combined baryonic and dark matter components of *Eris* produce a total gravitational potential shape of $q_z = 0.88$, and that this shape allows for the formation of tidal streams that show remarkable consistency with the Pal 5 and GD-1 data. Using Bayesian inference, we demonstrate that the *Eris* model works as well as the best fit (spherical) potential presented by Bovy et al. (2016). Remarkably, the tidal streams Pal 5 and GD-1 can be well-reproduced by the potential of *Eris* without any rescaling or modifications, suggesting that our current models of ΛCDM structure formation including dissipative physics can produce galaxies that have total halo potentials in substantial agreement with current observational constraints available for the Milky Way.

We find that the total potential *Eris* works comparably as well as the spherical NFW halo model found by Bovy et al. (2016), even though the dark matter halo in *Eris* has an axis ratio of $c/a = 0.65$. The baryonic dissipation alters the dark matter halo shape from prolate to an oblate halo aligned with the disk, consistent with previous studies (e.g., Kazantzidis et al. 2004; Gustafsson et al. 2006; Abadi et al. 2010). The more centrally concentrated mass distribution of *Eris* relative to an NFW halo allows for *Eris* to permit realistic stream orbits, despite its flattened shape. The *Eris* halo shape is substantially rounder than in the equivalent dark matter–only simulation *ErisDark*; we find that *ErisDark* cannot reproduce the observed stellar stream data.

With the ongoing release of the *Gaia* data, we expect exciting additional probes of the Milky Way potential that will further test our models for ΛCDM galaxy formation. *Gaia* will provide a full-mission proper motion precision of $\sigma_\mu \lesssim 100 \mu\text{as yr}^{-1}$ for
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Time evolution of the *Eris* dark matter halo. The left and middle panels show the minor–major and intermediate–major axis ratios of the *Eris* halo potential as a function of time, split by the dark matter (dark blue), baryons (light blue), and total potential (medium blue) components. The axis ratios are measured at *r* = 15 kpc, close to where Pal 5 and GD-1 are located. The right panel displays the evolution of mass within 20 kpc of the halo for the dark matter (dark blue) and baryons (light blue). These results suggest that the halo remains stable over the last ~4 Gyr.</figcaption></figure>

single stars down to 20th magnitude (de Bruijne 2012). With multiple available stars at magnitudes *m* < 20 in each stream, the constraints on both the Pal 5 and GD-1 proper motions should sharpen substantially beyond the $\sigma_{\mu} \approx 150$ µas yr⁻¹ uncertainties from the current data (Koposov et al. 2010; Fritz & Kallivayalil 2015). The ultimate precision afforded by *Gaia* should distinguish between models like *Eris* and more spherical halos, if the final stream proper motion constraints reach $\sigma_{\mu} \lesssim 50$ µas yr⁻¹.

We thank the anonymous referee for comments that helped improve our work. We also thank Annalisa Pillepich for providing access to the *ErisDark* simulation. Support for this work was provided by the Chinese Scholarship Council (B.D.), and by NASA through grant HST-AR-13904.001-A (P.M.) and contract NNG16PJ25C (B.E.R.). P.M. also thanks the Préfecture of the Ile-de-France Region for the award of a Blaise Pascal International Research Chair, managed by the Fondation de l'Ecole Normale Supérieure. This research utilized the NSF-funded Hyades supercomputer at UCSC.

## ORCID iDs

Biwei Dai https://orcid.org/0000-0003-2734-0294

Brant E. Robertson https://orcid.org/0000-0002-4271-0364

Piero Madau https://orcid.org/0000-0002-6336-3293

## References

* Abadi, M. G., Navarro, J. F., Fardal, M., Babul, A., & Steinmetz, M. 2010, MNRAS, 407, 435

* Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27

* Bovy, J. 2014, ApJ, 795, 95

* Bovy, J., Bahmanyar, A., Fritz, T. K., & Kallivayalil, N. 2016, ApJ, 833, 31

* Bovy, J., Bird, J. C., García Pérez, A. E., et al. 2015, ApJ, 800, 83

* Bovy, J., Nidever, D. L., Rix, H.-W., et al. 2014, ApJ, 790, 127

* Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115

* Bowden, A., Belokurov, V., & Evans, N. W. 2015, MNRAS, 449, 1391

* Bowden, A., Evans, N., & Williams, A. 2016, MNRAS, 460, 329

* Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40

* Brooks, A. M., & Zolotov, A. 2014, ApJ, 786, 87

* Buote, D. A., & Canizares, C. R. 1998, MNRAS, 298, 811

* Buote, D. A., Jeltema, T. E., Canizares, C. R., & Garmire, G. P. 2002, ApJ, 577, 183

* Chan, T., Kereš, D., Óñorbe, J., et al. 2015, MNRAS, 454, 2981

* Clutton-Brock, M. 1973, Ap&SS, 23, 55

* Cole, S., & Lacey, C. 1996, MNRAS, 281, 716

* Debattista, V. P., Moore, B., Quinn, T., et al. 2008, ApJ, 681, 1076 de Bruijne, J. H. J. 2012, Ap&SS, 341, 31

* Diemand, J., Kuhlen, M., Madau, P., et al. 2008, Natur, 454, 735

* Dotter, A., Sarajedini, A., & Anderson, J. 2011, ApJ, 738, 74

* Dubinski, J. 1994, ApJ, 431, 617

* Dubinski, J., & Carlberg, R. 1991, ApJ, 378, 496

* Fellhauer, M., Belokurov, V., Evans, N., et al. 2006, ApJ, 651, 167

* Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

* Flores, R. A., & Primack, J. R. 1994, ApJL, 427, L1

* Flynn, C., Sommer-Larsen, J., & Christensen, P. R. 1996, MNRAS, 281, 1027

* Frenk, C. S., White, S. D., Davis, M., & Efstathiou, G. 1988, ApJ, 327, 507

* Fritz, T. K., & Kallivayalil, N. 2015, ApJ, 811, 123

* Gnedin, O. Y., Kravtsov, A. V., Klypin, A. A., & Nagai, D. 2004, ApJ, 616, 16

* Governato, F., Zolotov, A., Pontzen, A., et al. 2012, MNRAS, 422, 1231

* Grillmair, C. J., & Dionatos, O. 2006, ApJL, 643, L17

* Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742, 76

* Guedes, J., Mayer, L., Carollo, M., & Madau, P. 2013, ApJ, 772, 36

* Gustafsson, M., Fairbairn, M., & Sommer-Larsen, J. 2006, PhRvD, 74, 123522

* Hammer, F., Puech, M., Chemin, L., Flores, H., & Lehnert, M. D. 2007, ApJ, 662, 322

* Harris, W. E. 1996, AJ, 112, 1487

* Helmi, A. 2004, ApJL, 610, L97

* Hernquist, L. 1990, ApJ, 356, 359

* Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375

* Hoekstra, H., Yee, H. K., & Gladders, M. D. 2004, ApJ, 606, 67

* Ibata, R., Lewis, G. F., Irwin, M., Totten, E., & Quinn, T. 2001, ApJ, 551, 2941

* Ibata, R. A., Lewis, G. F., Irwin, M. J., & Quinn, T. 2002, MNRAS, 332, 915

* Ibata, R. A., Lewis, G. F., & Martin, N. F. 2016, ApJ, 819, 1

* Ibata, R. A., Lewis, G. F., Thomas, G., Martin, N. F., & Chapman, S. 2017, ApJ, 842, 120

* Jing, Y., & Suto, Y. 2002, ApJ, 574, 538

* Johnston, K. V., Law, D. R., & Majewski, S. R. 2005, ApJ, 619, 800

* Kazantzidis, S., Abadi, M. G., & Navarro, J. F. 2010, ApJL, 720, L62

* Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., et al. 2004, ApJL, 611, L73
* King, I. 1962, AJ, 67, 471

* Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82

* Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260

* Kuhlen, M., Diemand, J., & Madau, P. 2007, ApJ, 671, 1135

* Küpper, A. H. W., Balbinot, E., Bonaca, A., et al. 2015, ApJ, 803, 80

* Küpper, A. H. W., Lane, R. R., & Heggie, D. C. 2012, MNRAS, 420, 2700

* Kuzma, P. B., Da Costa, G. S., Keller, S. C., & Maunder, E. 2015, MNRAS, 446, 3297

* Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229

* Ledoux, V., & Van Daele, M. 2016, ACM Transactions on Mathematical Software (TOMS), 42, 29

* Loebman, S. R., Ivezić, Ž., Quinn, T. R., et al. 2014, ApJ, 794, 151

* Madau, P., Shen, S., & Governato, F. 2014, ApJL, 789, L17

* Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533

* Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJL, 524, L19
Moore, B. 1994, *Natur*, **370**, 629

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, **462**, 563

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, **490**, 493

Ngan, W., Bozek, B., Carlberg, R. G., et al. 2015, ApJ, **803**, 75

Ngan, W., Carlberg, R. G., Bozek, B., et al. 2016, ApJ, **818**, 194

Odenkirchen, M., Grebel, E. K., Dehnen, W., et al. 2003, AJ, **126**, 2385

Odenkirchen, M., Grebel, E. K., Dehnen, W., Rix, H.-W., & Cudworth, K. M. 2002, AJ, **124**, 1497

Odenkirchen, M., Grebel, E. K., Rockosi, C. M., et al. 2001, ApJL, **548**, L165

Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, MNRAS, **454**, 2092

Parker, L. C., Hoekstra, H., Hudson, M. J., Van Waerbeke, L., & Mellier, Y. 2007, ApJ, **669**, 21

Pearson, S., Küpper, A. H. W., Johnston, K. V., & Price-Whelan, A. M. 2015, ApJ, **799**, 28

Pillepich, A., Kuhlen, M., Guedes, J., & Madau, P. 2014, ApJ, **784**, 161

Pontzen, A., & Governato, F. 2012, MNRAS, **421**, 3464

Portail, M., Wegg, C., Gerhard, O., & Martinez-Valpuesta, I. 2015, MNRAS, **448**, 713

Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Numerical Recipes. The Art of Scientific Computing (Cambridge: Cambridge Univ. Press)

Rahmati, A., & Jalali, M. A. 2009, MNRAS, **393**, 1459

Robin, A. C., Marshall, D. J., Schultheis, M., & Reylé, C. 2012, A&A, **538**, A106

Sandford, E., Küpper, A. H. W., Johnston, K. V., & Diemand, J. 2017, MNRAS, **470**, 522

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, **403**, 1829

Smith, M. C., Evans, N. W., & An, J. H. 2009, ApJ, **698**, 1110

Sokolowska, A., Babul, A., Mayer, L., Shen, S., & Madau, P. 2017, arXiv:1708.07820

Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, **373**, 1074

Tissera, P. B., White, S. D., Pedrosa, S., & Scannapieco, C. 2010, MNRAS, **406**, 922

Vera-Ciro, C. A., Sales, L. V., Helmi, A., et al. 2011, MNRAS, **416**, 1377

Vivas, A. K., & Zinn, R. 2006, AJ, **132**, 714

Wadsley, J. W., Stadel, J., & Quinn, T. 2004, NewA, **9**, 137

Warren, M. S., Quinn, P. J., Salmon, J. K., & Zurek, W. H. 1992, ApJ, **399**, 405

Weinberg, M. D. 1999, AJ, **117**, 629

Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJL, **827**, L23

Xue, X. X., Rix, H. W., Zhao, G., et al. 2008, ApJ, **684**, 1143

Zhao, H. 1996, MNRAS, **278**, 488