SPENDING TOO MUCH TIME AT THE GALACTIC BAR: CHAOTIC FANNING OF THE OPHIUCHUS STREAM

The Ophiuchus stream (Bernard et al. 2014; Sesar et al. 2015) is a recently discovered stellar tidal stream that sits above the Galactic bulge at a Galactocentric radius and height (R, z) ≈ (1.5, 4.3) kpc. All observational evidence suggests that the stream is a completely disrupted globular cluster: The stream stars have (1) a small positional dispersion orthogonal to the extended direction of the stream (width ≈10 pc, length ≈1.5 kpc); (2) no detectable over-density along the stream that could be the progenitor system; (3) a small velocity dispersion ≈0.4 km s⁻¹; and (4) an old stellar population (≈12 Gyr) estimated from isochrone fitting (Sesar et al. 2015, hereafter S15).

There are a number of peculiarities about the observed kinematics of the Ophiuchus stream. For example, the de-projected length of the visible part of the stream is short given the age of its stellar population (≈1.5 kpc). S15 fit an orbit to the kinematics of the stream stars in a static, axisymmetric model for the gravitational field of the inner Galaxy and ran N-body simulations of globular clusters on this orbit. S15 find that—on this orbit—the portion of the stream visible as an over-density in main-sequence stars must have been formed in the last ≲400 Myr for the stream to remain as short as it is observed. This dynamical age is at odds with the old (≈10–12 Gyr) stellar population: The abrupt end of the stream suggests that the cluster apparently fully disrupted at once in the last 400 Myr. Another puzzle is the existence of blue horizontal branch (BHB) stars close to the stream (within a few degrees) with similar radial velocities, but with a large dispersion in both sky position and velocity (Sesar et al. 2016, hereafter S16). The stream has a very distinct and large line-of-sight velocity (≈290 $\mathrm{km}\quad {{\rm{s}}}^{-1}$) and is therefore easily detected above the background halo population. Four BHB stars have been detected with line-of-sight velocities >230 $\mathrm{km}\quad {{\rm{s}}}^{-1}$ that lie close to an extrapolation of the stream on the sky. This makes them likely members of the stream, as their velocities are in stark contrast to the background halo population (see Section 4, S16). Yet, they have a velocity dispersion ≈75 times larger than the measured internal velocity dispersion of the stream stars. These stars hint at the existence of associated low-density, high-dispersion features that were not modeled in S15 and are not predicted by the N-body simulations from this prior work, which assume a sudden, total disruption of the stream progenitor.

The orbit fit and N-body simulations in S15 used a static, axisymmetric potential to represent the Milky Way potential, but it is well-known that the Galactic bulge contains a triaxial, rotating, bar-like structure several kpc in size (e.g., Blitz & Spergel 1991; Weinberg 1992; Dwek et al. 1995; Wegg & Gerhard 2013). Given the proximity of the stream to the center of the Galaxy, the time-dependent, triaxial potential of the Galactic bar must be taken into account when modeling the orbit of the Ophiuchus stream. The presence of a bar-like perturbation to the potential will change the orbit of the stream progenitor and the orbit structure in the inner galaxy (Zotos 2012; Gajda et al. 2015; Portail et al. 2015a). Bar-like features can also introduce a significant number of chaotic orbits in their vicinity (Weinberg 2015) and generate resonances that may also affect stream formation (Hattori et al. 2015).

Recent work has shown that dynamical chaos can dramatically alter the density evolution of tidal streams (e.g., Fardal et al. 2015; Price-Whelan et al. 2016). Along certain chaotic orbits, the stream stars will spread much faster in 3D position than from ordinary phase-mixing and, depending on the orbital phase at which the stream is observed, may develop large, low-density "fans" of stars at the ends of a stream (Pearson et al. 2015; Price-Whelan et al. 2016). As a first application of this theoretical understanding, we study whether stream-fanning—chaotic or simply from density evolution in a triaxial, time-dependent potential—could plausibly explain the observed properties of the Ophiuchus stream. In particular, we consider whether such models can reproduce:

• 1.  
  the apparent shortness and fast density truncation of the stream;
• 2.  
  the increased positional dispersion of the four new candidate members from S16;
• 3.  
  the large velocity dispersion of the S16 stars.

We do not aim to perfectly represent the observed data, but rather to explore the plausibility of explaining the peculiarities of the stream using chaotic stream fanning. Note that an alternate model was recently proposed that instead places the stream progenitor on an orbit in resonance with the bar (Hattori et al. 2015). We discuss the differences between these two models in Section 5.

In Section 2 we describe the methods used in this work: in Section 2.1 we describe the models we use for the gravitational potential of the Galaxy, in Section 2.2 we outline the probabilistic procedure we use to fit orbits to the stream data (explained in detail in Appendix B), and in Section 2.3 we explain the simple method we use to generate mock streams. In Section 3 we discuss the results from fitting orbits to the data in a static, axisymmetric potential model and several potential models with a time-dependent bar. In Section 4 we generate mock streams along these orbits to argue that chaotic stream-fanning is a plausible explanation for the observational peculiarities of the Ophiuchus stream. We discuss the implications of this work and possibilities for future work in Section 5 and we conclude in Section 6.

Our goal is to (1) assess whether the Galactic bar can produce chaotic orbits in the vicinity of the Ophiuchus stream and (2) determine if chaotic density evolution of tidal debris stripped from the progenitor of this stream can explain the apparent shortness of the stream and low-density, high-dispersion stars beyond the extent of main sequence stars observed in PS1. In this section, we describe the potential models we use to represent the galaxy and outline the methods we use to detect and quantify the strength of chaos for individual orbits. We then describe the likelihood function we use for fitting orbits to the stream stars. Finally, we describe how we generate mock stellar streams for a progenitor on a given orbit.

Throughout we assume the Sun is at Galactocentric position (x, y, z) = (−8.3, 0, 0) kpc (e.g., Schönrich 2012) with velocity $({v}_{x},{v}_{y},{v}_{z})=(-11.1,250,7.25)\;\mathrm{km}\quad {{\rm{s}}}^{-1}$ (e.g., Schönrich et al. 2010; Schönrich 2012).

2.1. Potential Models

To integrate orbits and to compute chaos indicators we must choose a gravitational potential model to represent the potential of the Milky Way. The key feature of the potential that we would like to capture is the time-dependence and triaxiality of the Galactic bar. Recent work has used stellar number counts of Red Clump giant stars in the Galactic bulge to constrain dynamical models of the bar (Portail et al. 2015b). Measurements of the total mass of the bar feature from this study are largely consistent with past work (e.g., Wang et al. 2012), however the measured pattern speed and present bar angle are significantly discrepant and this difference is not fully understood. We construct a parametrized potential model consisting of a triaxial, time-dependent (rotating) bulge component added to simple models for the disk and halo of the Milky Way. We describe below how we fix the parameters of the disk, halo, and bar or bulge component, but explore different choices for the time-dependence and orientation of the bar. We also define a static potential with a spherical bulge for comparison.

These potential models are meant to be representative rather than definitive. The uncertainty in the Milky Way potential within Galactocentric radii of r ≲ 4 kpc and outside of r ≳ 15 kpc are large enough that trying to match the exact density distribution of the Ophiuchus stream is not a useful exercise. Instead, we consider qualitatively different potentials that allow us to isolate and study the affect of chaotic stream-fanning of tidal debris in the vicinity of the stream.

2.1.1. Barred Potential

We use a spherical Navarro–Frenk–White potential to represent the dark matter halo (Navarro et al. 1996) parametrized as

$$\begin{eqnarray}&&{\rm{\Phi }}(r)=-{v}_{h}^{2}\;\displaystyle \frac{\mathrm{ln}(1+r/{r}_{s})}{r/{r}_{s}}\end{eqnarray} \tag{ 1 }$$

and a Miyamoto–Nagai potential for the disk (Miyamoto & Nagai 1975). For the bar component, we use a basis function expansion (BFE) of the potential and density of the bar with expansion coefficients derived for a triaxial, exponential bar density (Wang et al. 2012, hereafter W12). We use the pre-computed expansion coefficients used in W12, which were computed from a low-order expansion of the triaxial bar density used in Dwek et al. (1995).³ We have implemented the BFE computation of the potential, density, and gradient of the potential in C and Python and the code is publicly available on GitHub.⁴

The BFE representation fixes the axis ratios of the bar—that is, the exponential scale lengths along the three axes of the bar were adopted from Dwek et al. (1995) when the expansion coefficients were calculated in W12; all other potential parameter values are given in Table 1. The mass of the halo is fixed and the mass of the disk and bar are varied in order to qualitatively reproduce the flatness and amplitude of the circular velocity curve of the Milky Way (Bovy et al. 2012). Figure 1, top left shows the circular velocity along the line connecting the Sun to the Galactic center in this model (the Galactic x axis). Figure 1, top right shows contours of constant surface density for a face-on (left) and edge-on (right) view of this potential model with the bar angle set to 20° (compare to, e.g., Figure 3 in Portail et al. 2015b). We consider a grid of nine parameter combinations of bar angle and pattern speed. Model names and parameter values are given in Table 2.

Zoom In Zoom Out Reset image size

Table 1.  The Disk Potential Scale Lengths (a, b) were Adopted following (Bovy 2015) to Match the Exponential Scale Length of the Disk (Bovy & Rix 2013) and Local Dark-matter Density (e.g., Bovy & Tremaine 2012)

Component	Parameter	Value
Disk	Mdisk	4 × 1010 ${M}_{\odot }$
	a	3 kpc
	b	0.28 kpc
Spheroid	Msph	5 × 109 ${M}_{\odot }$
	c	0.2
Halo	vc	185.8 $\mathrm{km}\quad {{\rm{s}}}^{-1}$
	rs	30 kpc
Bar	Mbar	1.8 × 1010 ${M}_{\odot }$

Note. The halo mass scale is set by specifying the circular velocity at the scale radius, v_c, and the scale velocity in Equation (1) is given by ${v}_{h}^{2}={v}_{c}^{2}/(\mathrm{ln}2-1/2)$. The bar mass is taken from recent 3D density modeling of red clump stars in the galactic bulge (Portail et al. 2015b). The other bar parameters are listed in Table 2 next to the corresponding model name.

Download table as:  ASCIITypeset image

Table 2.  Present-Day Bar Angle (α) and Pattern Speed (Ω_p) for the Nine Parameter Pairs Considered in this Work

Name	α (deg)	Ωp (km s−1 kpc−1)
bar1	20	40
bar2	20	50
bar3	20	60
bar4	25	40
bar5	25	50
bar6	25	60
bar7	30	40
bar8	30	50
bar9	30	60

Note. These values span the range of recent measurements from a variety of techniques (Dwek et al. 1995; Wang et al. 2012, 2013; Wegg & Gerhard 2013).

Download table as:  ASCIITypeset image

2.1.2. Static Potential

For comparison, we also define a time-independent potential model with a purely spherical bulge. In this model, we set the bar mass to 0 and instead add a spheroidal component represented with a Hernquist potential (Hernquist 1990). Parameters for this potential model are given in Table 3. Figure 1, bottom left shows the circular velocity along the line connecting the Sun to the Galactic center in this model (the Galactic x axis). Figure 1, bottom right shows contours of constant surface density for a face-on (left) and edge-on (right) view of this potential model.

Table 3.  Same as Table 1, Except: the Disk Mass is Increased to Account for Removing the Bar Component, a Spheroidal Bulge Component is Added

Component	Parameter	Value
Disk	Mdisk	6 × 1010 ${M}_{\odot }$
	a	3 kpc
	b	0.28 kpc
Halo	vc	185.8 $\mathrm{km}\quad {{\rm{s}}}^{-1}$
	rs	30 kpc
Spheroid	Msph	1.2 × 1010 ${M}_{\odot }$
	c	0.3

Download table as:  ASCIITypeset image

2.2. Fitting Orbits to the Ophiuchus Stream

2.3. Generating Mock Streams

Figures 2–3 summarize our results from fitting orbits to the BHB stream stars. Shown in each panel are the high-probability Ophiuchus stream stars (black points, to which orbits are fit) and orbits integrated from samples from the posterior probability over orbital parameters (blue lines). The four "fanned" BHB stars from S16 are shown as gray squares and are not included in the orbit fitting procedure. We only show one of the barred potentials: the end-to-end integration time of the orbit over the observed extent of the stream is only ≈6 Myr, so the derived orbits are extremely similar in all potentials (the time-dependence of the bar potential is not significant over such short timescales). As was previously noted, there is a marginal discrepancy between the orbit fits and the observed proper motions, but it is possible there is a systematic offset in the proper motion measurements (see Figure 10 in Sesar et al. 2015). We conclude that it is too soon to tell whether there is in fact an inconsistency. The orbital periods are typically ≈170–200 Myr with pericenters r_p ≈ 4 kpc and apocenters r_a ≈ 12–15 kpc. Though the coordinate and velocity parameter values in observed coordinates are very similar between each potential model, the resulting orbits are quite different. For the posterior samples in each potential, we take the mean values of the coordinate and velocity parameters at ϕ₁ = 0 and convert to Galactocentric coordinates (e.g., Table 4). Figure 4 shows projections of these "mean" orbits in each potential model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Estimated Mean and Standard Deviation of Samples from the Marginal Posterior Distributions over Each Parameter in Our Orbit Fit Model At ϕ₁ = 0 (the Posterior Distributions are Very Close to Gaussian)

Name	ϕ2 (deg)	d⊙ (kpc)	μl (mas yr−1)	μb (mas yr−1)	vr (km s−1)
static	−0.03 ± 0.05	8.3 ± 0.05	−7.4 ± 0.1	0.9 ± 0.1	288.9 ± 0.9
bar1–9	−0.03 ± 0.05	8.35 ± 0.05	−7.4 ± 0.1	0.9 ± 0.1	289.0 ± 1.0

sϕ2 (deg)	sd (kpc)	${s}_{{v}_{r}}$ (km s−1)
0.20 ± 0.04	0.31 ± 0.10	2.9 ± 0.8
0.21 ± 0.05	0.31 ± 0.14	3.2 ± 0.9

Note. For the barred potentials, all mean values are the same because the time-dependence of the bar does not impact the orbit fit over the short length of the stream. We have made samples from the full posterior distribution available with this article and provide code to transform to and from stream coordinates (see http://adrian.pw/ophiuchus for more information).

Download table as:  ASCIITypeset image

For the posterior samples in each potential, we also compute the maximum Lyapunov exponent (${\lambda }_{{\rm{max}}}$) and corresponding Lyapunov time (${t}_{\lambda }=1/{\lambda }_{{\rm{max}}}$) to assess whether each orbit is chaotic. For strongly chaotic orbits, the Lyapunov time is still an appropriate indicator of chaos and of the timescale over which chaos is important for tidal debris (Price-Whelan et al. 2016). Figure 5 shows distributions of Lyapunov times for orbits drawn from the posterior distributions from orbit fitting in each potential model. All orbits in the static potential have Lyapunov times ${t}_{\lambda }\gt 20\;{\rm{Gyr}}$ and we consider them to be regular (no panel is shown for these orbits). All orbits sampled from each barred potential are strongly chaotic with Lyapunov times that range from ${t}_{\lambda }\approx 400$–1100 Myr. From the left column to the right column, the distribution of Lyapunov times shifts slightly towards lower values suggesting that the orbits around Ophiuchus are more strongly chaotic for larger pattern speeds (within the range considered in this work).

Zoom In Zoom Out Reset image size

The orbits sampled from the orbit-fit posteriors in the barred potentials are all strongly chaotic. We have tried computing the frequency diffusion rate for these orbits as an independent check of their chaotic timescale but have found that, over consecutive integration windows, the frequency recovery fails or is unreliable because the frequency spectrum changes dramatically over timescales of ≈10–20 orbital periods. To understand whether time-dependence or the triaxiality of the bar is the dominant cause of chaos, we have also performed the orbit-fit and Lyapunov time experiments in a potential with a fixed bar (Ω_p = 0) at the same bar angles explored above. We find that the Lyapunov times for orbits fit to the stream in these potentials are consistent with being regular. Computing the frequency diffusion rate instead shows that these orbits are likely weakly chaotic: the frequency diffusion times (Price-Whelan et al. 2016) are ≈10⁴ Gyr. We therefore conclude that the time-dependence of the bar is the primary source of chaos in the inner Galaxy.

Here we study whether the observed abrupt drop in density and possible fanned debris stars can be explained without assuming a sudden change in the mass-loss history of the cluster. In particular, we are interested in whether the mock streams formed around the strongly chaotic progenitor orbits in the barred potentials can explain these features while having been steadily disrupted over many Gigayears.

For each potential model, we randomly sample 256 orbits from the orbital parameter posterior distributions and generate mock streams along each orbit. We use the method outlined above (Section 2.3) to generate the streams and set the free parameters as follows: (1) we evolve the progenitors for 6 Gyr along each orbit prior to the current position to explore stream models where the shortness of the stream is not due to an instantaneous disruption 400 Myr ago, (2) we release star particles every 0.5 Myr (uniformly in time) to densely sample the final density distribution, (3) we set the progenitor mass to m = 10⁴ ${M}_{\odot }$(as was estimated by S15), and (4) we set the disruption time of the progenitor equal to the last time at which a pericenter and disk crossing coincide (in each case, this is at t ≈ −200 Myr). After the disruption time, we continue releasing the star particles uniformly in time (every 0.5 Myr) rather than releasing a "burst" of particles at once. We therefore expect that the density of the most recently disrupted debris will be systematically higher for the model streams as compared to the observed stream.

For each generated mock stream, we compute the likelihood of the data (now including all BHB stars from S15 and S16, e.g., all points in Figure 2) given the star particles by estimating the phase-space model density using a kernel density estimate with a Gaussian kernel (see, e.g., Bonaca et al. 2014). For each i data point ${{\boldsymbol{x}}}_{i}$ and each k model point with coordinates ${{\boldsymbol{y}}}_{k}$, we compute the likelihood by converting the model point position and velocity into heliocentric coordinates and evaluate

$$\begin{eqnarray}&&p({{\boldsymbol{x}}}_{i}\;| \;{{\boldsymbol{y}}}_{k},{{\boldsymbol{\sigma }}}_{i},{\boldsymbol{h}})={ \mathcal N }({{\boldsymbol{x}}}_{i}\;| \;{{\boldsymbol{y}}}_{k},{{\boldsymbol{\sigma }}}_{i}^{2}+{{\boldsymbol{h}}}^{2})\end{eqnarray} \tag{ 2 }$$

where ${ \mathcal N }(x\;| \;\mu ,{\sigma }^{2})$ is the normal distribution with mean μ and variance σ² and ${\boldsymbol{h}}$ represents the diagonal of the bandwidth matrix, ${\boldsymbol{H}}$, used for the density estimate, ${\boldsymbol{h}}={\rm{diag}}({\boldsymbol{H}})$. We fix the bandwidth parameters as follows

$$\begin{eqnarray}&&{\boldsymbol{h}}=\left(\begin{array}{c}0.2\;{\rm{deg}}\\ 0.2\;{\rm{deg}}\\ 0.2\;{\rm{kpc}}\\ 0.2\;\mathrm{km}\quad {{\rm{s}}}^{-1}\end{array}\right)\end{eqnarray} \tag{ 3 }$$

for sky position in Galactic coordinates (l, b), distance, and radial velocity. The full likelihood for all N data points given K model points is

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \prod _{i}^{N}\displaystyle \frac{1}{K}\displaystyle \sum _{k}^{K}\;p({{\boldsymbol{x}}}_{i}\;| \;{{\boldsymbol{y}}}_{k},{{\boldsymbol{\sigma }}}_{i},{\boldsymbol{h}}).\end{eqnarray} \tag{ 4 }$$

Figures 6–7 show the final particle positions and line-of-sight velocities in heliocentric coordinates for the maximum likelihood mock streams (points) in each potential. Vertical lines show the approximate extent of the part of the stream visible in main-sequence stars (excluding the BHB stars from S16). Figures 8–9 show the same for proper motion components. There are a few interesting features to note from these panels:

• 1.  
  Even in the static potential (Figure 6, leftmost column), there is a slight decrease in the density for the model stream towards higher Galactic longitudes, l. This is a projection effect: The portion of the stream at larger l is closer and points almost directly towards the Sun so that the debris covers a larger area on the sky.
• 2.  
  The density of the mock stream in the static potential decreases slowly rather than abruptly as is observed. This is more easily seen in Figure 10 where each contour level represents a factor of 10 difference in projected surface density. In the static potential model the length of dense debris extends much farther than the observed extent of the stream (vertical dashed lines), whereas in some of the barred potential models the stars released earlier have "fanned" and are associated with much lower density debris (e.g., bar8).
• 3.  
  The four high-dispersion BHB stars beyond the end of the stream (from S16) do not match in position and velocity with the particle distribution from even the maximum likelihood stream model in the static potential. In some barred potentials the chaotic evolution of the stream stars can lead to over-densities of stars with an increased positional dispersion and significantly discrepant velocities (e.g., bar8).
• 4.  
  None of the stream models—static or barred—produce an appreciable density of stars with line-of-sight velocities near the S16 BHB star with the largest velocity (≈320 $\mathrm{km}\quad {{\rm{s}}}^{-1}$). This star is either an interesting Ophiuchus member star or is associated with some other kinematic substructure.
• 5.  
  For the barred potentials, the stream morphology is very sensitive to the properties of the bar (especially the pattern speed) and to the initial conditions of the orbit. We have found that the morphology can vary significantly between nearby orbits in the same potential model (because these are strongly chaotic orbits), but the overall characteristics remain similar: along more strongly chaotic orbits, the debris "fans" more and the apparent dense part of the model streams is shorter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The density truncation of the mock streams in each potential model is more clearly seen in terms of the density contrast between stream stars and background stars, visualized in Figure 11: This figure shows mock sky-density maps of stars generated by superimposing the maximum likelihood model streams over a noisy background of stars. The number of mock stream star particles used to generate the map has been normalized such that the total number of stars within the observed extent of the stream (between 585 > l > 381) is equal to the number of stars attributed to the Ophiuchus stream in the PS1 data (N ≈ 500 Bernard et al. 2014).

Zoom In Zoom Out Reset image size

The viewing angle and stream geometry is more clearly demonstrated in Figure 12, which shows x–z projections of the star particles in Galactocentric Cartesian coordinates (gray) along with the position of the Sun (symbol at (x, z) = (−8.3, 0) kpc) and the "window" of the heliocentric, sky-position plots of Figures 6–7 (shown as blue lines).

Zoom In Zoom Out Reset image size

The model streams presented here do not reproduce all observed features of the Ophiuchus stream (the details of the potential and the orbit predict vastly different phase-space morphologies for the fanned part). Instead, these results illustrate that chaotic evolution of tidal debris can plausibly explain the peculiar features of the stream. If the cluster progenitor was on a regular orbit fit to the observed stream stars in a static, axisymmetric potential model for the Milky Way, it would have to have disrupted entirely within the last ≈300–600 Myr in order to explain the shortness and density profile. In addition, the four most recently identified BHB stars with similar distances and line-of-sight velocities would have to be (highly unlikely) chance alignments of halo stars. If instead the progenitor were on a chaotic orbit (because of the influence of the Galactic bar): (1) stars stripped early will have "fanned" out and would thus be harder to observe and (2) the nearby, high-velocity-dispersion BHB stars can be naturally explained by this chaotic stream-fanning. We consider the second scenario to be more plausible: our understanding of the formation and evolution of the Ophiuchus stream is that the progenitor object has been orbiting and steadily losing stars over the last several Gigayears, but only the stars stripped from the most recent few pericentric passages remain coherent enough to be detected as a stream-like over-density in the PS1 data.

It is still too early to say for sure that chaotic stream-fanning is occurring for the Ophiuchus stream. Deeper follow-up imaging and spectroscopy over a larger area region around the stream will be needed to test the predictions of this work and compare with other possible scenarios. For example, recent work has shown that if the Ophichus stream progenitor orbit is in resonance with the bar, the debris can remain short for at least 1 Gyr (Hattori et al. 2015). In their model, there would be no nearby, high-dispersion debris, and the pattern speed of the bar would be related to the orbital frequencies of the progenitor orbit. However, it has not been demonstrated whether this proposed scheme can explain the shortness of the stream over timescales closer to the age of the stellar population (≈10 Gyr). With more information about the density distribution of stars in the stream and better proper motion measurements we would be able to (1) help distinguish models for the Milky Way bar independent from current methods and (2) begin to model the survivability of globular clusters orbiting in the central Milky Way (e.g., Gnedin & Ostriker 1997).

5.1. Future Work: Modeling the Galactic Bar

5.2. Future Work: The Orbits and Survivability of Inner Milky Way Globular Clusters

We have shown that, with a qualitative but observationally motivated potential model for the Galactic bar, the orbits of the Ophiuchus stream stars are likely sensitive to the time-dependence and shape of the bar potential. For modeling the stream density itself, it is therefore crucial to include this component of the Galactic potential. By fitting orbits to kinematic data for members of the stream in Milky Way-like potential models, we have found that orbits in the vicinity of the Ophiuchus stream are strongly chaotic for a range of bar parameters (pattern speeds and present-day angles). Using mock stellar stream models generated assuming a globular cluster-mass progenitor object, we have shown that the apparent shortness of the stream and the existence of nearby stars with very high velocity dispersion are plausibly explained by chaotic density evolution of the stars stripped from the progenitor object.

This is the first time chaos has been used to explain the morphology of a stellar stream and the first observational evidence for the importance of chaos in the Galactic halo. It also highlights the importance of including the Galactic bar in dynamical modeling of the Milky Way's inner halo and has important implications for future modeling of streams near in this region. With more Ophiuchus stream members, density and velocity information over a larger region near the stream, and better models for the internal structure of the Galactic bar, careful modeling of this stream could lead to tight constraints on the structure and evolution of the bar.

APW is supported by a National Science Foundation Graduate Research Fellowship under Grant No. 11-44155. KVJ and APW acknowledge support from NSF grant AST-1312196. This material is based on work supported by the National Aeronautics and Space Administration under Grant No. NNX15AK78G issued through the Astrophysics Theory Program. APW acknowledges the staff at the MPIA for their support and assistance. The authors wish to acknowledge Victor Debattista, Melissa Ness, and Sarah Pearson for useful discussions. APW acknowledges David Bowie (1947–2016) for his continuous inspiration and artistic vision. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013). This work additionally relied on Columbia University's Yeti compute cluster, and we acknowledge the Columbia HPC support staff for assistance.

The transformation matrix is approximately represented as

$$\begin{eqnarray*}&&{\left(\begin{array}{c}x\\ y\\ z\end{array}\right)}_{{\rm{Oph}}}\\ \approx \;\left(\begin{array}{ccc}0.84922096554 & 0.07001279040 & 0.52337554476\\ -0.27043653641 & -0.79364259852 & 0.54497294023\\ 0.45352820359 & -0.60434231606 & -0.65504391727\end{array}\right)\\ &&\quad \times \;{\left(\begin{array}{c}x\\ y\\ z\end{array}\right)}_{{\rm{Gal}}}\end{eqnarray*}$$

but the precise transformation and coordinate frame is implemented in Python using the Astropy coordinates package.⁵ This code is hosted on GitHub.⁶

Our goal is to infer the posterior probability distributions over orbital initial conditions, ${{\boldsymbol{w}}}_{0}={(l,b,{\rm{DM}},{\mu }_{l},{\mu }_{b},{v}_{r})}_{0}$, given a potential, Φ, and kinematic data for each i stream star, ${{\boldsymbol{x}}}_{i}={(l,b,{\rm{DM}},{\mu }_{l},{\mu }_{b},{v}_{r})}_{i}$. In this notation, (l, b) are Galactic coordinates, ${\rm{DM}}$ is the distance modulus, $({\mu }_{l},{\mu }_{b})$ are proper motions in the Galactic frame, and v_r is the radial velocity. We assume that the sky coordinates for each star are known perfectly well (have zero uncertainty) and transform the data to a rotated, heliocentric coordinate system that is aligned with the stream and centered on the median sky position of the BHB stars in the densest part of the stream (all BHB stars except the "fanned" stars: cand15, cand26, cand49, cand54 from Sesar et al. 2016). We represent the longitude and latitude in these coordinates as (ϕ₁, ϕ₂) and the rotation matrix to transform from Galactic to these coordinates is given in Appendix A. We treat the stream longitude, ϕ₁, as the perfectly known, independent variable so that all other coordinates can be expressed as functions of this longitude (e.g., ϕ₂(ϕ₁), DM(ϕ₁), etc.). This methodology is similar to that used in Koposov et al. (2010) and Sesar et al. (2015).

B.1. Likelihood

We include three nuisance parameters in our likelihood to account for the internal dispersion of the stream: in observed coordinates, these are the on-sky positional dispersion, ${s}_{{\phi }_{2}}$, a distance (modulus) dispersion, s_DM, and a radial velocity dispersion, ${s}_{{v}_{r}}$ (the proper motion uncertainties are sufficiently large that we can not resolve the velocity dispersion in these coordinates).⁷ We add two additional nuisance parameters for controlling the amount of time to integrate forwards, t_f, and backwards, t_b, from the given initial conditions, which ultimately controls the length of the section of orbit that is compared to the stream star data. For brevity in the equations below, we define ${\boldsymbol{s}}=({s}_{{\phi }_{2}},{s}_{{\rm{DM}}},{s}_{{v}_{r}})$, ${{\boldsymbol{\sigma }}}_{i}={({\sigma }_{{\rm{DM}}},{\sigma }_{{\mu }_{l}},{\sigma }_{{\mu }_{b}},{\sigma }_{{v}_{r}})}_{i}$, and ${\boldsymbol{\theta }}=({{\boldsymbol{w}}}_{0},{\rm{\Phi }},{t}_{b},{t}_{f})$.

For a given set of initial conditions (${{\boldsymbol{w}}}_{0}$), we compute a model orbit as follows: (1) transform the initial conditions to Galactocentric coordinates, (2) integrate the orbit forward and backward by t_f and t_b, respectively, in the potential Φ, (3) transform all orbit points (time-steps) back to observed coordinates, and (4) define interpolating functions for each coordinate as a function of stream longitude, ϕ₁, using cubic splines—e.g., functions ${\tilde{\phi }}_{2}({\phi }_{1})$, $\tilde{{\rm{DM}}}({\phi }_{1})$, $\tilde{{\mu }_{l}}({\phi }_{1})$, $\tilde{{\mu }_{b}}({\phi }_{1})$, $\tilde{{v}_{r}}({\phi }_{1})$. These functions let us compute the predicted values of each of these coordinates at the longitudes of each observed star, ϕ_{1, i}.

We assume that each observed kinematic component is independent so that the likelihood of the data for a given star, ${{\boldsymbol{x}}}_{i}$, with uncertainties, ${{\boldsymbol{\sigma }}}_{i}$, is given by the product over the likelihoods for each dimension of the data:

$$\begin{eqnarray}p({{\boldsymbol{x}}}_{i}\;| \;{{\boldsymbol{\sigma }}}_{i},{\boldsymbol{s}},{\boldsymbol{\theta }}) & = & p({\phi }_{2,i}\;| \;{\phi }_{1,i},{s}_{{\phi }_{2}},{\boldsymbol{\theta }})\\ & & \times p({{\rm{DM}}}_{i}\;| \;{\phi }_{1,i},{\sigma }_{{\rm{DM}},i},{s}_{{\rm{DM}}},{\boldsymbol{\theta }})\\ & & \times p({\mu }_{l,i}\;| \;{\phi }_{1,i},{\sigma }_{{\mu }_{l},i},{\boldsymbol{\theta }})\\ & & \times p({\mu }_{b,i}\;| \;{\phi }_{1,i},{\sigma }_{{\mu }_{b},i},{\boldsymbol{\theta }})\\ & & \times p({v}_{r,i}\;| \;{\phi }_{1,i},{\sigma }_{{v}_{r},i},{s}_{{v}_{r}},{\boldsymbol{\theta }}).\end{eqnarray} \tag{ 5 }$$

The uncertainties in these observed coordinate components are assumed to be normally distributed away from the model values: using the notation

$$\begin{eqnarray}&&{ \mathcal N }(x\;| \;\mu ,{\sigma }^{2})=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^{2}}}\;\mathrm{exp}\left(-\displaystyle \frac{{(x-\mu )}^{2}}{2{\sigma }^{2}}\right)\end{eqnarray} \tag{ 6 }$$

the likelihoods are

$$\begin{eqnarray}&&p({\phi }_{2,i}\;| \;{\phi }_{1,i},{s}_{{\phi }_{2}},{\boldsymbol{\theta }})={ \mathcal N }({\phi }_{2,i}\;| \;{\tilde{\phi }}_{2}({\phi }_{1,i}),{s}_{{\phi }_{2}}^{2})\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&p({{\rm{DM}}}_{i}\;| \;{\phi }_{1,i},{\sigma }_{{\rm{DM}},i},{s}_{{\rm{DM}}},{\boldsymbol{\theta }})\\ &&\quad ={ \mathcal N }({{\rm{DM}}}_{i}\;| \;\tilde{{\rm{DM}}}({\phi }_{1,i}),{s}_{{\rm{DM}}}^{2}+{\sigma }_{{\rm{DM}},i}^{2})\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&p({\mu }_{l,i}\;| \;{\phi }_{1,i},{\sigma }_{{\mu }_{l},i},{\boldsymbol{\theta }})={ \mathcal N }({\mu }_{l,i}\;| \;\tilde{{\mu }_{l}}({\phi }_{1,i}),{\sigma }_{{\mu }_{l,i}}^{2})\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&p({\mu }_{b,i}\;| \;{\phi }_{1,i},{\sigma }_{{\mu }_{b},i},{\boldsymbol{\theta }})={ \mathcal N }({\mu }_{b,i}\;| \;\tilde{{\mu }_{b}}({\phi }_{1,i}),{\sigma }_{{\mu }_{b,i}}^{2})\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&p({v}_{r,i}\;| \;{\phi }_{1,i},{\sigma }_{{v}_{r},i},{s}_{{v}_{r}},{\boldsymbol{\theta }})={ \mathcal N }({v}_{r,i}\;| \;\tilde{{v}_{r}}({\phi }_{1,i}),{s}_{{v}_{r}}^{2}+{\sigma }_{{v}_{r},i}^{2}).\end{eqnarray} \tag{ 11 }$$

We assume the data from each star is independent and identically distributed so that the full likelihood is the product over the likelihoods for all N stars:

$$\begin{eqnarray}&&p(\{{{\boldsymbol{x}}}_{i}\}\;| \;\{{{\boldsymbol{\sigma }}}_{i}\},{\boldsymbol{s}},{\boldsymbol{\theta }})=\displaystyle \prod _{i}^{N}p({{\boldsymbol{x}}}_{i}\;| \;{{\boldsymbol{\sigma }}}_{i},{\boldsymbol{s}},{\boldsymbol{\theta }}).\end{eqnarray} \tag{ 12 }$$

B.2. Priors

For the intrinsic dispersion parameters, we use logarithmic (scale-invariant) priors such that p(s) ∝ s⁻¹. For the integration time parameters, we use uniform priors, ${ \mathcal U }(a,b)$ (over the range a–b),

$$\begin{eqnarray}&&p({t}_{f})={ \mathcal U }(1,100)\;{\rm{Myr}}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&p({t}_{b})={ \mathcal U }(-100,-1)\;{\rm{Myr}}.\end{eqnarray} \tag{ 14 }$$

Note that present-day is t = 0. For computational efficiency, we place strong priors on the minimum and maximum longitudes of the model points, (ϕ_1,min, ϕ_1,max) so that the model orbit does not integrate for longer than necessary. In particular, we set

$$\begin{eqnarray}&&p({\phi }_{1,{\rm{min}}}\;| \;{\boldsymbol{\theta }})={ \mathcal N }({\phi }_{1,{\rm{min}}}\;| \;\mathrm{min}({\phi }_{1,i}),{s}_{{\phi }_{2}}^{2})\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&p({\phi }_{1,{\rm{max}}}\;| \;{\boldsymbol{\theta }})={ \mathcal N }({\phi }_{1,{\rm{max}}}\;| \;\mathrm{max}({\phi }_{1,i}),{s}_{{\phi }_{2}}^{2}).\end{eqnarray} \tag{ 16 }$$

For the orbital initial condition components, we use uniform priors in each cartesian position component over the range (−200, 200) kpc. For velocity, we use a Gaussian prior on the magnitude of the total velocity, v, with a dispersion of 150 $\mathrm{km}\quad {{\rm{s}}}^{-1}$,

$$\begin{eqnarray}&&{ \mathcal N }(v\;| \;0,{(150\mathrm{km}{{\rm{s}}}^{-1})}^{2}).\end{eqnarray} \tag{ 17 }$$

We keep the potential, Φ, fixed. In total, this model has 10 parameters (5 phase-space coordinates, 5 nuisance parameters).

The full expression for the posterior probability, $p({\boldsymbol{s}},{{\boldsymbol{w}}}_{0},{t}_{b},{t}_{f}\;| \;\{{{\boldsymbol{x}}}_{i}\},\{{{\boldsymbol{\sigma }}}_{i}\},{\rm{\Phi }})$, is the joint likelihood (Equation (12)) multiplied by all priors described above (Equations (13)–(17)).