Spending too much time at the Galactic bar: chaotic fanning of the Ophiuchus stream
Abstract
The Ophiuchus stellar stream is peculiar: (1) its length is short given the age of its constituent stars, and (2) several probable member stars have dispersions in sky position and velocity that far exceed those seen within the stream. The stream’s proximity to the Galactic center suggests that its dynamical history is significantly influenced by the Galactic bar. We explore this hypothesis with models of stream formation along orbits consistent with Ophiuchus’ properties in a Milky Way potential model that includes a rotating bar. In all choices for the rotation parameters of the bar, orbits fit to the stream are strongly chaotic. Mock streams generated along these orbits qualitatively match the observed properties of the stream: because of chaos, stars stripped early generally form lowdensity, highdispersion “fans” leaving only the most recently disrupted material detectable as a strong overdensity. Our models predict that there should be a significant amount of lowsurfacebrightness tidal debris around the stream with a complex phasespace morphology. The existence of or lack of these features could provide interesting constraints on the Milky Way bar and would rule out formation scenarios for the stream. This is the first time that chaos has been used to explain the properties of a stellar stream and is the first demonstration of the dynamical importance of chaos in the Galactic halo. The existence of long, thin streams around the Milky Way, presumably formed along non or weaklychaotic orbits, may represent only a subset of the total population of disrupted satellites.
1 Introduction
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 . 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 overdensity along the stream that could be the progenitor system; (3) a small velocity dispersion 0.4 ; 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 deprojected 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 Nbody simulations of globular clusters on this orbit. S15 find that—on this orbit—the portion of the stream visible as an overdensity in mainsequence 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 Myrs. 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 lineofsight velocity (290 ) and is therefore easily detected above the background halo population. Four BHB stars have been detected with lineofsight velocities 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 lowdensity, highdispersion features that were not modeled in S15 and are not predicted by the body simulations from this prior work, which assume a sudden, total disruption of the stream progenitor.
The orbit fit and body simulations in S15 used a static, axisymmetric potential to represent the Milky Way potential, but it is wellknown that the Galactic bulge contains a triaxial, rotating, barlike 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 timedependent, triaxial potential of the Galactic bar must be taken into account when modeling the orbit of the Ophiuchus stream. The presence of a barlike perturbation to the potential will change the orbit of the stream progenitor and the orbit structure in the inner galaxy (Zotos, 2012; Portail et al., 2015a; Gajda et al., 2015). Barlike 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; PriceWhelan et al., 2016). Along certain chaotic orbits, the stream stars will spread much faster in 3D position than from ordinary phasemixing and, depending on the orbital phase at which the stream is observed, may develop large, lowdensity “fans” of stars at the ends of a stream (Pearson et al., 2015; PriceWhelan et al., 2016). As a first application of this theoretical understanding, we study whether streamfanning—chaotic or simply from density evolution in a triaxial, timedependent potential—could plausibly explain the observed properties of the Ophiuchus stream. In particular, we consider whether such models can reproduce:

the apparent shortness and fast density truncation of the stream;

the increased positional dispersion of the four new candidate members from S16;

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 timedependent bar. In Section 4 we generate mock streams along these orbits to argue that chaotic streamfanning 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.
2 Methods
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 lowdensity, highdispersion 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 (e.g., Schönrich, 2012) with velocity (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 timedependence 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, timedependent (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 timedependence 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 and outside of 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 streamfanning of tidal debris in the vicinity of the stream.
Barred potential
We use a spherical NavarroFrenkWhite potential to represent the dark matter halo (Navarro et al., 1996) parametrized as
(1) 
and a MiyamotoNagai 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 precomputed expansion coefficients used in W12, which were computed from a loworder expansion of the triaxial bar density used in Dwek et al. (1995).
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 axis). Figure 1, top right shows contours of constant surface density for a faceon (left) and edgeon (right) view of this potential model with the bar angle set to (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.
Component  Parameter  Value 

Disk  
3 kpc  
0.28 kpc  
Spheroid  
0.2  
Halo  185.8  
30 kpc  
Bar 
Name  [deg]  [] 

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 
Static potential
For comparison, we also define a timeindependent 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 axis). Figure 1, bottom right shows contours of constant surface density for a faceon (left) and edgeon (right) view of this potential model.
Component  Parameter  Value 

Disk  
3 kpc  
0.28 kpc  
Halo  185.8  
30 kpc  
Spheroid  
0.3 
2.2 Fitting orbits to the Ophiuchus stream
In each of the potentials described above, we fit orbits to the measured kinematics of BHB stars that are highlikelihood members of the Ophiuchus stream (Sesar et al., 2015, 2016). The details of this procedure and a definition of the likelihood function we use are presented in Appendix B. We use an ensemble Markov Chain Monte Carlo (MCMC) algorithm (Goodman & Weare, 2010) implemented in Python (emcee) to generate samples from the posterior distribution over the parameters in our orbitfitting model (ForemanMackey et al., 2013). The algorithm uses an ensemble of individual “walkers” to adapt to the geometry of the parameterspace being explored. In all cases, we use 80 walkers (8 times the number of parameters).
To initialize these walkers, we first run an optimization routine to maximize the likelihood: We use the Powell algorithm implemented in Scipy (Powell, 1964; Jones et al., 2001–) to minimize the negative, loglikelihood. To generate initial conditions for the walkers, we sample from Gaussian distributions centered on the maximum likelihood values. For the coordinates, we set the dispersions of these Gaussians to 1/1000 of the median uncertainties of the stars. For the nuisance parameters, we set the dispersions to 1/1000 of their maximum likelihood values.
For each potential, we run the MCMC walkers for a burnin period of 512 steps and then reinitialize the walkers from their positions at the end of this run. This erases any relics of the initialization procedure outlined above. After burnin, we run the walkers for an additional 512 steps. For each parameter, we compute the autocorrelation times, , of the Markov chains and thin the chains by taking every sample. This reduces the number of samples, but ensures that our posterior samples are effectively independent.
2.3 Generating mock streams
To generate mock stellar streams, we use a method similar to that presented in Fardal et al. (2015): Star particles are ‘released’ from a progenitor system near the Lagrange points with a dispersion in position and velocity that is set by the mass and orbit of the progenitor. We draw samples from the posterior probability distributions over orbital parameters from fitting orbits to the stream star members and use these as the progenitor orbital parameters (Section B). For a given progenitor orbit—the 6D position of the orbit today—we integrate the orbit backwards in time for a given integration period. From the endpoint of the backwardsintegration (e.g., the past position), we begin integrating the orbit forward in time, but now at each timestep a star particle is released near each of the Lagrange points of the progenitor. This simplification assumes that the massloss rate is constant in time, which is not assumed in Fardal et al. (2015) but has been shown to closely match body simulations of globular cluster disruption (Küpper et al., 2012). The position of the Lagrange points and the scale of the dispersion in position and velocity are set by the progenitor mass, . The star particles are drawn from Gaussians centered on the Lagrange points (in position) and the progenitor (in velocity) and the full parametrization of the release distribution is given in Fardal et al. (2015). This method has been shown to reproduce the morphologies of body simulations of stellar streams, but requires far less computing time because it relies only on integrating testparticle orbits.
We make one modification to this method based on the idea that the Ophiuchus stream progenitor has been fullydisrupted. We add an additional parameter to the stream generation routine to specify the time of disruption, . At this time, we set the offset of the Lagrange point to 0: In terms of the parametrization in (Fardal et al., 2015), we set and but preserve the dispersion in the release radius and velocity. That is, any star released after is drawn from a Gaussian with positional and velocity dispersion fixed to their respective values at with no offset from the progenitor orbit. This massloss history is intended to mimic the expected gradual evaporation of a globular cluster over a tidallylimited boundary (i.e. driven by twobody relaxation and gravitational shocks over many Gigayears) with final disruption likely to occur once the tidal boundary is less than the core radius of the cluster. The physics of this disruption is not followed exactly but rather the disruption rate and final disruption time are set by the hypothesis that the most recent combined pericenter and disk shock fully disrupted the cluster, but, critically, that the cluster has been losing debris over its entire orbital history.
3 Results I: Orbit fits and chaos
Figures 2–3 summarize our results from fitting orbits to the BHB stream stars. Shown in each panel are the highprobability 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 grey squares and are not included in the orbit fitting procedure. We only show one of the barred potentials: The endtoend 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 timedependence 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 and apocenters –. 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 and convert to Galactocentric coordinates (e.g., Table 4). Figure 4 shows projections of these “mean” orbits in each potential model.
For the posterior samples in each potential, we also compute the maximum Lyapunov exponent (MLE, ) and corresponding Lyapunov time () 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 (PriceWhelan 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 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 –. 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).
The orbits sampled from the orbitfit 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 timedependence or the triaxiality of the bar is the dominant cause of chaos, we have also performed the orbitfit and Lyapunov time experiments in a potential with a fixed bar () 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 (PriceWhelan et al., 2016) are . We therefore conclude that the timedependence of the bar is the primary source of chaos in the inner Galaxy.
4 Results 2: Stream models for the Ophiuchus stream
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 massloss 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 Myrs 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 (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 ). 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.
name  [deg]  [kpc]  [mas yr]  [mas yr]  [km s] 

static  
bar1–9 
[deg]  [kpc]  [km s] 

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 phasespace model density using a kernel density estimate with a Gaussian kernel (see, e.g., Bonaca et al., 2014). For each data point and each model point with coordinates , we compute the likelihood by converting the model point position and velocity into heliocentric coordinates and evaluate
(2) 
where is the normal distribution with mean and variance and represents the diagonal of the bandwidth matrix, , used for the density estimate, . We fix the bandwidth parameters as follows
(3) 
for sky position in Galactic coordinates (, ), distance, and radial velocity. The full likelihood for all data points given model points is
(4) 
Figures 6–7 show the final particle positions and lineofsight 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 mainsequence stars (excluding the BHB stars from S16). Figure 8–9 show the same for proper motion components. There are a few interesting features to note from these panels:

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, . This is a projection effect: The portion of the stream at larger is closer and points almost directly towards the Sun so that the debris covers a larger area on the sky.

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).

The four highdispersion BHB stars beyond the end of the stream (from S16) don’t 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 overdensities of stars with an increased positional dispersion and significantly discrepant velocities (e.g., bar8).

None of the stream models—static or barred—produce an appreciable density of stars with lineofsight velocities near the S16 BHB star with the largest velocity (). This star is either an interesting Ophiuchus member star or is associated with some other kinematic substructure.

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.
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 skydensity 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 ) is equal to the number of stars attributed to the Ophiuchus stream in the PS1 data ( Bernard et al., 2014).
The viewing angle and stream geometry is more clearly demonstrated in Figure 12, which shows  projections of the star particles in Galactocentric Cartesian coordinates (grey) along with the position of the Sun (symbol at ) and the “window” of the heliocentric, skyposition plots of Figures 6–7 (shown as blue lines).
5 Discussion
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 phasespace 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 lineofsight 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, highvelocitydispersion BHB stars can be naturally explained by this chaotic streamfanning. 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 streamlike overdensity in the PS1 data.
It is still too early to say for sure that chaotic streamfanning is occurring for the Ophiuchus stream. Deeper followup 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, highdispersion 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
The current kinematic data for the Ophiuchus BHB stars suggests that the stream is sensitive to the gravitational potential of the Milky Way bar—with better measurements of the velocities and a larger sample of member stars we will constrain parameters for the bar model. Current measurements of the pattern speed, angle, and structure of the Milky Way’s bulge and bar are largely discrepant (e.g., Wang et al., 2012, 2013; Wegg & Gerhard, 2013; Antoja et al., 2014). Most of the methods that infer these quantities rely on modeling the density or kinematics of stars at low Galactic latitudes and must therefore handle challenges with completeness and dust extinction. The Ophiuchus stream offers a unique opportunity to independently measure these quantities by modeling the density and kinematics of stars associated with the stream.
5.2 Future work: The orbits and survivability of inner Milky Way globular clusters
If the Ophiuchus stream formed from a globular cluster on a chaotic orbit and we happen to be witnessing its final demise, what does this imply about the population of clusters that have already been fully disrupted? The existence of strongly chaotic orbits in this region would limit the expected number of cold stellar streams in the inner Galaxy and enhance the rate of mixing of the debris. The fraction of strongly chaotic orbits that would lead to chaotic fanning or fast dispersal of tidal debris should therefore be related to the amount of kinematic substructure in the inner Galaxy. Indeed, first suggestions of kinematic substructure in the bulge have been found, but further modeling is needed to understand whether these hints could be signature of a widely dispersed globular cluster population. If so, a stronger theoretical understanding of the prevalence of these features could be combined with future kinematic surveys (from, e.g., Gaia) to place constraints on longstanding puzzles about the primordial globular cluster population (e.g., Murali & Weinberg, 1997; Gnedin & Ostriker, 1997).
6 Conclusions
We have shown that, with a qualitative but observationallymotivated potential model for the Galactic bar, the orbits of the Ophiuchus stream stars are likely sensitive to the timedependence 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 Waylike 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 presentday angles). Using mock stellar stream models generated assuming a globular clustermass 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.
Appendix A Transformation from Galactic to Ophiuchus stream coordinates
The transformation matrix is approximately represented as
but the precise transformation and coordinate frame is implemented in Python using the Astropy coordinates package.
Appendix B Fitting orbits to stellar streams
Our goal is to infer the posterior probability distributions over orbital initial conditions, , given a potential, , and kinematic data for each stream star, . In this notation, are Galactic coordinates, is the distance modulus, are proper motions in the Galactic frame, and 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 perfectlyknown, independent variable so that all other coordinates can be expressed as functions of this longitude (e.g., , , 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 onsky positional dispersion, , a distance (modulus) dispersion, , and a radial velocity dispersion, (the proper motion uncertainties are sufficiently large that we can’t resolve the velocity dispersion in these coordinates).
For a given set of initial conditions (), we compute a model orbit as follows: (1) transform the initial conditions to Galactocentric coordinates, (2) integrate the orbit forward and backward by and , respectively, in the potential , (3) transform all orbit points (timesteps) back to observed coordinates, and (4) define interpolating functions for each coordinate as a function of stream longitude, , using cubic splines—e.g., functions , , , , . These functions let us compute the predicted values of each of these coordinates at the longitudes of each observed star, .
We assume that each observed kinematic component is independent so that the likelihood of the data for a given star, , with uncertainties, , is given by the product over the likelihoods for each dimension of the data:
(B1) 
The uncertainties in these observed coordinate components are assumed to be normally distributed away from the model values: using the notation
(B2) 
the likelihoods are
(B3)  
(B4)  
(B5)  
(B6)  
(B7) 
We assume the data from each star is independent and identically distributed (i.i.d.) so that the full likelihood is the product over the likelihoods for all stars:
(B8) 
b.2 Priors
For the intrinsic dispersion parameters, we use logarithmic (scaleinvariant) priors such that . For the integration time parameters, we use uniform priors, (over the range –),
(B9)  
(B10) 
Note that presentday is . For computational efficiency, we place strong priors on the minimum and maximum longitudes of the model points, so that the model orbit does not integrate for longer than necessary. In particular, we set
(B11)  
(B12) 
For the orbital initial condition components, we use uniform priors in each cartesian position component over the range . For velocity, we use a Gaussian prior on the magnitude of the total velocity, , with a dispersion of ,
(B13) 
We keep the potential, , fixed. In total, this model has 10 parameters (5 phasespace coordinates, 5 nuisance parameters).
Footnotes
 affiliation: Department of Astronomy, Columbia University, 550 W 120th St., New York, NY 10027, USA
 affiliation: To whom correspondence should be addressed: adrn@astro.columbia.edu
 affiliation: MaxPlanckInstitut für Astronomie, Königstuhl 17, D69117 Heidelberg, Germany
 affiliation: Department of Astronomy, Columbia University, 550 W 120th St., New York, NY 10027, USA
 affiliation: MaxPlanckInstitut für Astronomie, Königstuhl 17, D69117 Heidelberg, Germany
 The coefficients presented in W12 are for just the cosine terms (the in Hernquist & Ostriker (1992) or the in Lowing et al. (2011)) because all sine terms have zero coefficients for a triaxial density function.
 https://github.com/adrn/biff
 See http://adrian.pw/ophiuchus for more information
 https://github.com/adrn/ophiuchus
 We assume that the dispersion in these coordinates is constant over the observed (short) section of the stream. This may be a bad assumption.
References
 Antoja, T., Helmi, A., Dehnen, W., et al. 2014, A&A, 563, A60
 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
 Bernard, E. J., Ferguson, A. M. N., Schlafly, E. F., et al. 2014, MNRAS, 443, L84
 Blitz, L., & Spergel, D. N. 1991, ApJ, 379, 631
 Bonaca, A., Geha, M., Küpper, A. H. W., et al. 2014, ApJ, 795, 94
 Bovy, J. 2015, ApJS, 216, 29
 Bovy, J., & Rix, H.W. 2013, ApJ, 779, 115
 Bovy, J., & Tremaine, S. 2012, ApJ, 756, 89
 Bovy, J., Allende Prieto, C., Beers, T. C., et al. 2012, ApJ, 759, 131
 Dwek, E., Arendt, R. G., Hauser, M. G., et al. 1995, ApJ, 445, 716
 Fardal, M. A., Huang, S., & Weinberg, M. D. 2015, MNRAS, 452, 301
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
 Gajda, G., Lokas, E. L., & Athanassoula, E. 2015, ArXiv eprints, arXiv:1511.04253
 Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223
 Goodman, J., & Weare, J. 2010, Comm. App. Math. Comp. Sci., 5, 65
 Hattori, K., Erkal, D., & Sanders, J. L. 2015, ArXiv eprints, arXiv:1512.04536
 Hernquist, L. 1990, ApJ, 356, 359
 Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375
 Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python, [Online; accessed 20150211]
 Koposov, S. E., Rix, H.W., & Hogg, D. W. 2010, ApJ, 712, 260
 Küpper, A. H. W., Lane, R. R., & Heggie, D. C. 2012, MNRAS, 420, 2700
 Lowing, B., Jenkins, A., Eke, V., & Frenk, C. 2011, MNRAS, 416, 2697
 Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
 Murali, C., & Weinberg, M. D. 1997, MNRAS, 291, 717
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
 Pearson, S., Küpper, A. H. W., Johnston, K. V., & PriceWhelan, A. M. 2015, ApJ, 799, 28
 Portail, M., Wegg, C., & Gerhard, O. 2015a, MNRAS, 450, L66
 Portail, M., Wegg, C., Gerhard, O., & MartinezValpuesta, I. 2015b, MNRAS, 448, 713
 Powell, M. J. D. 1964, CompJ, 7, 155
 PriceWhelan, A. M., Johnston, K. V., Valluri, M., et al. 2016, MNRAS, 455, 1079
 Schönrich, R. 2012, MNRAS, 427, 274
 Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829
 Sesar, B., Bovy, J., Bernard, E. J., et al. 2015, ApJ, 809, 59
 Sesar, B., PriceWhelan, A. M., Cohen, J. G., et al. 2016, ApJ, 816, L4
 Wang, Y., Mao, S., Long, R. J., & Shen, J. 2013, MNRAS, 435, 3437
 Wang, Y., Zhao, H., Mao, S., & Rich, R. M. 2012, MNRAS, 427, 1429
 Wegg, C., & Gerhard, O. 2013, MNRAS, 435, 1874
 Weinberg, M. D. 1992, ApJ, 384, 81
 —. 2015, ArXiv eprints, arXiv:1508.05959
 Zotos, E. E. 2012, RAA, 12, 500