{'fraction_non_alphanumeric': 0.02940944881889764, 'fraction_numerical': 0.012992125984251968, 'mean_word_length': 4.6534609392388155, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Probabilistic model for constraining the Galactic potential using tidal streams

We present a generative probabilistic model for a tidal stream and demonstrate how this model is used to constrain the Galactic potential. The model takes advantage of the simple structure of a stream in angle and frequency space for the correct potential. We investigate how the method performs on full 6D mock stream data, and mock data with outliers included. As currently formulated the technique is computationally costly when applied to data with large observational errors, but we describe several modifications that promise to make the technique computationally tractable.


INTRODUCTION
A key goal in the study of the Milky Way is mapping the dark matter distribution of the Galaxy. Locally this is achieved through dynamical measurements of the disc components (Kuijken & Gilmore 1989;Garbari et al. 2012;Bovy & Tremaine 2012;Zhang et al. 2013), whilst tidal streams present a tantalising prospect for constraining the dark matter distribution on a more global scale. These long, filamentary structures are the remnants of satellites tidally disrupted by the Milky Way. The stars in the stream are now essentially orbiting freely in the Galactic potential, but importantly they form a dynamically coherent structure, such that their properties are closely linked together. By inspecting the phase-space structure of tidal streams it should be possible to infer properties of the Galactic potential, and, in particular, the dark matter distribution.
Many methods have been proposed to achieve this. The simplest, and therefore most fruitful to-date, are orbit-fitting techniques (Binney 2008;Eyre & Binney 2009a,b;Koposov et al. 2010;Deg & Widrow 2014). Several authors have acknowledged the problems with simply fitting a stream with a single orbit (Eyre & Binney 2011;Sanders & Binney 2013a;Lux et al. 2013). Therefore, much work has been done on constructing methods for recovering the potential without assuming the stream delineates an orbit. One group of methods seeks the potential which minimises the spread in the integrals of motion (Peñarrubia et al. 2012). Others allow for the expected range of orbits present in a stream given the progenitor properties (Johnston et al. 1999;Varghese et al. 2011;Price-Whelan & Johnston 2013).
The dynamics of tidal streams are very simply expressed in angle-action coordinates (Tremaine 1999 Eyre & Binney 2011) denoted as (J , θ) throughout the paper. In Sanders & Binney (2013b) we demonstrated the power of using the frequencies and angles of the stream to construct a measure of the goodness of fit of the potential. The frequencies, Ω, are the derivatives of the Hamiltonian, H, with respect to the actions i.e. Ω = ∂H/∂J . The best potential was deemed the one in which the angle and frequency distributions are aligned. This method suffered several drawbacks: we were doing the inference in model space, not the data space; it was difficult to assess the errors in the obtained potential parameters; it was awkward to handle errors in the data, and the method did not behave well for large errors.
Data for tidal streams have large associated errors and are expected to be contaminated with many stars that are not stream members. Streams lie at large distances from the Sun, far out in the halo. As such, the associated errors in the distance can be significant, and small proper motion errors translate into large transverse velocity errors. For many streams we don't have full 6D data. Additionally, a stream is identified as a filament in the observable space -often l and b. Stream members are then extracted by cutting appropriately in these coordinates, and cleaning up by using the additional observables. This process obviously introduces many outliers to the stream data, which are perhaps just members of a smoother background halo, whilst also potentially throwing out stars which are members of the stream.
This state of affairs leads us to analyse the data by constructing probabilistic models that correctly handle large errors, missing data, and contaminants. Such an approach is much more robust than previous efforts, and lends itself perfectly to being combined with other independent measurements of the Galactic potential. Here we present a probabilistic model for tidal streams that may be used to infer the properties of the Galactic potential. The model is expressed in the space of observables, but relies heavily on the ex-pected structure of streams in angle-action space. In Section 2 we motivate our choice of model by considering an idealised case of a Gaussian structure in angle-action space evolving in time. In Section 3 we use these insights to write down a practical model for the stream, and discuss how it may be used to infer properties of the Galactic potential. In Section 4 we infer the parameters of a simple two-parameter potential from mock stream observations using our model.

FORMALISM & MOTIVATION
Given a set of observations, D, of N stars believed to be members of a stream, what can we infer about the Galactic potential? For star i we have a 6D set of Galactic coordinates Li = (l, b, s, v || , µ) with associated gaussian error distributions described by the standard deviations σi. Note that we can fit any missing data into this formalism by taking the associated error with the data point to be infinite.
Given the data we want to know the posterior distribution of the potential, Φ, given by p(Φ|D). Using Bayes theorem we have where p(Φ) is the prior on the potential, and the evidence p(D) is not important for the present exercise. We wish to evaluate the likelihood p(D|Φ). The probability of the data given the potential is related to the properties of the stream progenitor, C . This contains information about the current phase-space coordinates of the progenitor (i.e. the progenitor's actions, J 0, frequencies, Ω0, and angles, θ0), as well as the size (and internal properties) of the progenitor. Therefore, we write p(D|Φ) = dC p(C )p(D|Φ, C ), where the Jacobian factor is given by and the (x, v) coordinates are related to the Galactic coordinates in the usual way. We want to work with actions, angles and frequencies. Throughout the paper we use the Stäckel-fitting algorithm presented in Sanders (2012) to estimate these quantities. This algorithm fits a Stäckel potential to the region of the potential an orbit probes, and estimates the actions, angles and frequencies in the true potential as those in the best-fitting Stäckel potential. Therefore, we write where the angles and frequencies are related to (x, v) via the potential Φ using the Stäckel-fitting approximation and the Jacobian is given by where we have used the fact that (J , θ) are canonical coordinates, such that the phase-space volume is conserved under the transformation, and introduced the Hessian matrix D defined as This matrix can be calculated analytically in a Stäckel potential, so we extend the Stäckel-fitting algorithm of Sanders (2012) to estimate D. We give details of this in the appendix. We proceed by splitting p(θi, Ωi|Φ, C ) into two components To proceed further we must consider what we know about stream formation in angle-action space (Tremaine 1999). For each star in the stream we have where ti is the time since the particle was stripped from the progenitor and ∆θi(0) is the separation between the i th particle and the progenitor when the particle is released. D0 is the Hessian from equation (7) evaluated at the progenitor actions, J 0.
To motivate our choice of model we begin by assuming that J i is distributed as an isotropic Gaussian such that where a gives the spread of the action-space distribution. This is related to the mass: a ∝ M 2/3 . Such a simple model for the actionspace distribution is unrealistic (Eyre & Binney 2011) but our understanding of this simplistic model will aid in the construction of a more realistic model. Similarly we assume that ∆θi(0) is distributed as an isotropic Gaussian such that This Gaussian model for a stream in actions and initial angles was studied by Helmi & White (1999). From equation (9) the frequency is linearly related to the actions via the Hessian, D0, so we can write down the distribution for the frequencies as where as the spread in actions is small we have approximated the Jacobian by its value at the progenitor actions. This distribution is a multivariate normal distribution with principal axes along the principal eigenvectors of D0 and with width given by the corresponding eigenvalues. D0 is a symmetric matrix so has real eigenvalues and orthogonal eigenvectors. Note here that for long thin streams to form D0 has one eigenvalue much greater than the other two. In Sanders & Binney (2013a) we demonstrated that in a realistic Galactic potential this condition was satisfied for a large volume of action space. Therefore we write where λj andêj are the eigenvalues and eigenvectors of D −1 0 , and we have λ1 λ2 > λ3. Therefore we find that The distribution of the frequencies is 2D Gaussian perpendicular to a straight line in frequency space defined byê1. In the simple model presented here the distribution alongê1 is also Gaussian. However, we will later adopt a superior distribution alongê1 which better reflects the stream distribution. Now we must address the angle distribution. The angles depend upon the additional variables, ti and ∆θi(0). Therefore, we write (0)).

(15)
Given a time since stripping, a frequency separation, and an initial angle separation, the present angle separation is completely determined by equation (9) so Substituting this and equation (11) into equation (15) and performing the integral over ∆θi(0) using the δ-function we have that Now we rearrange the argument of the exponential such that and note that ∆Ωi ≈ λ1ê1(ê1 · ∆J i) so and Therefore, equation (17) becomes The first part is a 2D Gaussian perpendicular to the eigenvectorê1 (as with the frequencies) whilst the second part depends on when the particles were stripped from the progenitor and only affects the angle distribution along the vectorê1 i.e. (∆θi ·ê1). If we assume that p(ti) is uniform and ∆θ ∆θ(0) (i.e. the particle was stripped long enough ago that the time part of equation (9) dominates the initial angle separation from the progenitor -this is an assumption we made in Sanders & Binney (2013b)), we may perform the ti integral to find where the condition ensures that the stripping time for each stream member is positive, and less than some maximum stripping time, tmax. This expression demonstrates explicitly that the distribution perpendicular toê1 is independent of the distribution alongê1. In conclusion in this model both the angle and frequency space distributions are highly elongated along the vectorê1. This validates the procedure followed in Sanders & Binney (2013b).

MODEL
In the formalism of the previous section we made several assumptions that, while useful for illustrative purposes, we would like to relax. The assumption of isotropic ∆J distribution is not valid as evidenced in Eyre & Binney (2011). For a general action-space distribution we would still expect a highly anisotropic frequency space distribution but the principal eigenvector of this distribution will not be that of the Hessian matrix D, but some other vectorn, with the vectorsd1 andd2 perpendicular to this. The intricacies of the action-space distribution will be reflected in the frequency distribution along the vectorn. Additionally the angle-space distribution will also be highly elongated along this directionn. Analogous to a combination of equation (14) and (21) we write where the functions Ki define the stream distribution along the vectorn in the angle or frequency space. The quantities sj and vj are the widths perpendicular to this vector in angle and frequency space, respectively. ωj and γj are related to the present frequency and angle coordinates of the progenitor. We define the angles φ and ψ such that n = (sin φ cos ψ, cos φ cos ψ, sin ψ), we choosê We now specify the functions Ki defining the stream distribution along the vectorn. In frequency space the distribution alongn consists of two separated peaks corresponding to the leading and trailing tails of the stream (see next section). For simplicity we assume that each of these peaks is Gaussian. The distribution in angle-space depends upon both the distribution in frequency space and the distribution of stripping times. As in equation (22), we make the simple first-order assumption of a uniform stripping time distribution such that the distribution in angle space along the stream given a frequency separation is also uniform between 0 and some maximum stripping time, tmax. Therefore, we write 2Ωs gives the separation between the Gaussian peaks alongn in frequency space. When equation (23) is combined with equation (2) and (5) we have completely specified our model. Given a set of 6D stream data with associated errors we can assess the likelihood of a given potential by evaluating the integral of equation (2). It is defined by 15 progenitor parameters given by C = {φ, ψ, γj, ωj, sj, vj, tmax, Ωs} and N potential parameters.

MCMC
We sample from the posterior using Markov chain Monte Carlo. We use an affine-invariant sampler implemented in the emcee package from Foreman-Mackey et al. (2013). For each of the following tests we use a group of 144 walkers, and vary the nuisance parameters C as well as the potential parameters. For all scale parameters (i.e. sj, vj) we use a logarithmic flat prior, whilst for the other parameters we use uniform flat priors. To perform the integral over the errors in the calculation of the likelihood we use the Vegas Monte Carlo integration algorithm (Lepage 1978) implemented in the Gnu Science Library. Our stream model is very narrow whilst the error distribution for each observable coordinate can be very broad. Therefore there is a very small region of the 4D integrand which has any support. Using an adaptive integration scheme such as Vegas means we can rapidly focus on this small region.

TESTS
We test the above procedure using particles taken from a stream simulation. The potential we work with is the two-parameter (N = 2) logarithmic potential given by The two parameters of this potential are Vc and q. We set Vc = 220 km s −1 and q = 0.9 for the simulation. We place a cluster at the apocentre of the orbit shown in Fig. 1. This orbit was chosen due to its similarity to the orbit of the GD-1 stream (Koposov et al. 2010). 11 th pericentre passage of the progenitor) using the code GYRFAL-CON (Dehnen 2000(Dehnen , 2002. We take the resulting distribution of particles, remove the progenitor remnant, and rotate the coordinate frame such that the Sun is placed at the same azimuthal angle as the progenitor. From the resulting particles we randomly select 30 particles that lie in the range −200 • < l < −140 • . The chosen sample of particles are shown in Fig. 1 and Fig. 2. For the tests shown below we include the observational errors by scattering the observed coordinates by the appropriate gaussian errors. In Fig. 3 we plot the angles and frequencies in the correct potential, along with a cross-section through the model specified in equation (23). The parameters of the model were chosen as those that maximised the likelihood. We see that in the correct potential the angle and frequency structure of the stream takes on a simple linear distribution. In Fig. 4 we plot the projections of the frequency histograms along the vectorsn,d1 andd2 for all particles in the stream from the simulation. We see that the distribution perpen- dicular ton is approximately Gaussian. The distribution alongn consists of two peaks consisting of the leading and trailing tail. Each peak is skewed such that there is a longer tail towards larger values of |∆Ω ·n|. The structure of these peaks was discussed by Johnston (1998). They can be considered as superpositions of Gaussians, one for each pericentric passage, which have decreasing separation from the progenitor and increasing amplitude as the cluster dissolves. By only sampling a portion of the stream track on the sky we have limited the range of available angles of the stream particles, but we expect that the range of frequencies sampled is fair. We will miss some high frequency separations which are stripped earliest. In this situation the tmax parameter tells us about the time since the first of the observed particles were stripped.

No contaminants
We test the above method by considering a dataset where each of our stream particles has full 6D data, (l, b, s, v || , µ). We first demonstrate that the method works for error-free data. In Figure 5 we show the posterior distributions of the potential parameters. The correct potential is recovered with ∼ 0.5 km s −1 and 0.005 errors in Vc and q respectively.
Next, we consider an optimistic dataset where we assume each of the stars is an RR Lyrae observed by Spitzer (Price-Whelan & Johnston 2013). We adopt 2 per cent distance errors (σs/s), 5 km s −1 line-of-sight velocity errors (σ || ), and 0.1 mas yr −1 proper motion errors (σµ). Such a dataset is unrealistic for such a low mass stream. However, it is suitable for demonstrating the method. The posterior distribution of the recovered potential parameters is shown in Figure 6. We recover the correct potential parameters as Vc = (218.9 ± 1.3) km s −1 and q = (0.908 ± 0.015) approximating the posterior as an uncorrelated Gaussian.

Inclusion of outliers
Stream data are inevitably contaminated with stars that are not associated with the stream. Many authors attempt to remove these outliers by performing cuts in the observable space. However, it is much better to model the outliers. To include m outliers in our test we randomly select m stars from our simulation sample. We use the l and b values of these m stars for our outlier stars. For each star we sample a distance, s, from a uniform distribution between 6 kpc and 12 kpc. We convert the tuple (l, b, s) to a Cartesian position, x. At this position x we draw a set of velocities v from an isothermal distribution function given by with dispersion σ h = 100km s −1 .
In our model we assume that the likelihood of a star at angle- frequency coordinates (Ωi, θi) is given by where pS is the likelihood given it is a member of the stream as outlined above, and p h is the likelihood given that it is not a member of the stream (either a member of a smooth background population or another structure). is the probability of being an outlier, which is given a logarithmic uniform prior. We choose to specify the outlier model in angle-frequency space as it is simpler to normalize in these coordinates, and we are less sensitive to systematics arising from the Stäckel approximation. We set p h to be uniform in both the angle and frequency space. This model is clearly simplistic, and we are not correctly accounting for selection effects in the angles due to observing stars in some region of the Galaxy. To our 30-particle stream dataset we add 20 outliers from our halo model, such that true = 0.4, and consider the simple case where we have an error-free dataset. The input data-set is shown in Fig. 7 and the resulting posterior distributions for the potential parameters are given in Fig. 8. In Fig. 9 we show the posterior distribution for the outlier fraction. It peaks nicely around the input outlier fraction, but has fairly large scatter around this value. This is probably due to the simplistic nature of the background model, but also, as we can see in Fig. 7, there is significant overlap in observable space between the stream and outlier distributions. Further information such as metallicities would potentially provide a simpler way to disentangle the stream from the background.

DISCUSSION
The formalism presented in this paper provides the first steps towards constraining the potential of the Galaxy from a generative stream model expressed in angle-frequency space. Through the course of developing the necessary machinery it has become apparent that many improvements can be made. These are as follows: (i) The model for a cold stream is very narrow, whilst the error distribution for the observational quantities is expected to be large. The convolution of the error and model distributions only has support over a very narrow range in each observable coordinate. Here we have performed the integration using Vegas to concentrate rapidly on this small region. However, a more efficient scheme would be full forward modelling. To do this we would perform the Monte Carlo integration by sampling 'true' coordinates from the stream distribution and sum the appropriate Gaussian error distributions for each of the observed coordinates. For this we need an efficient scheme for going from a frequency-angle sample to observable space. The torus machine (McMillan & Binney 2008) was developed to convert from angle-actions to (x, v) efficiently. Such machinery seems ideal for our purpose. A torus is constructed by finding the coefficients of the generating function from a toy torus. We can construct a series of nearby tori which cover the small range in action space occupied by the stream particles. For a given set of angles we can find a corresponding (x, v) on each torus. For any action point which lies between tori we can construct an appropriate torus on-the-fly by interpolating the generating function coefficients of the nearby tori. A small complication here is that we wish to work with frequencies instead of actions so a scheme for finding J (Ω) is required. We expect that for the region of interest this function is simple and may be obtained by interpolation between tori. With this function we may also calculate the Hessian as a function of Ω.
(iii) The distributions in frequency and angle-space, Ki, along the principal stream direction,n, were taken in this paper to be simple Gaussians and a uniform distribution. This procedure is adequate for our purposes but more realistic distributions are required to reproduce the expected feathering in the stream (Küpper et al. 2012;Sanders & Binney 2013b). For instance we could adjust the distribution over stripping times to be more concentrated around pericentric passage of the progenitor.
(iv) Our choice of model parametrisation is perhaps more flexible than it should be. We have chosen to uncouple the widths of the Gaussians to give the model more freedom. In theory we could relate v1 to v2 via the Hessian matrix as in equation (14), and set s1 = s2 as in equation (21). Additionally all the widths are related to M 1/3 so there could be a way to link all four parameters into a single mass parameter. Also, we have allowed the stream to be oriented along some random direction,n. The direction depends on the exact structure of the action-space distribution, and the structure of the Hessian matrix. The first of these can be sensitive to internal structure in the progenitor. We gain information aboutn through the use of the angle and frequency structure as discussed in Sanders & Binney (2013b). However, we could instead choose to maken a function of the progenitor actions in the chosen potential. This would constrain the models in a more physically motivated fashion, and would provide more information when the quality of the data is reduced.
(v) We have analysed a simulated stream evolved in a nearspherical logarithmic potential. As discussed in Sanders & Binney (2013a) this does not exhibit substantial offset between the stream and orbit tracks. Such an offset is more apparent in flattened potentials. Therefore, more tests are required with more realistic potentials with disc components to validate the modelling approach presented here.
Recently Bovy (2014) presented a machinery very similar to that shown here for constructing models of tidal streams. He exploits the narrow range of frequencies in the stream to construct a simple linear map between (Ω, θ) and (x, v). The Gaussian structure in (Ω, θ)-space can then be simply translated into a Gaussian structure in (x, v) making marginalization over missing data simpler. Helmi & White (1999) present a similar approach to analysing a stream model consisting of a Gaussian structure in action-space and initial angle-space.

CONCLUSIONS
We have presented a probabilistic model for a tidal stream and used this to constrain the potential from a simulation. The stream model builds on the work of Sanders & Binney (2013b) by constructing a simple model for the stream in frequency and angle space. The presented formalism naturally accounts for the errors in stream data, and can also incorporate the possibility of stream data being contaminated with stars from a smooth halo population or another tidal structure. We have successfully recovered the potential parameters used to run an N-body simulation of a GD-1-like stream from errorfree data, data with small errors included, and data with outliers included.
As currently formulated the computational cost of implementing our approach increases significantly with the magnitude of the observational errors. We have described modifications that promise to mitigate this effect, and thus to make the approach a powerful technique for constraining the Galaxy's gravitational potential.