Properties of Dark Subhaloes from Gaps in Tidal Streams
Cold or Warm, the Dark Matter substructure spectrum must extend to objects with masses as low as , according to the most recent Lyman- measurements. Around a Milky Way-like galaxy, more than a thousand of these subhaloes will not be able to form stars but are dense enough to survive even deep down in the potential well of their host. There, within the stellar halo, these dark pellets will bombard tidal streams as they travel around the Galaxy, causing small but recognizable damage to the stream density distribution. The detection and characterization of these stream ruptures will allow us to constrain the details of the subhalo-stream interaction. In this work, for the first time, we will demonstrate how the properties of a subhalo, most importantly its mass and size, can be reliably inferred from the gap it produces in a tidal stream. For a range of realistic observational setups, mimicking e.g. SDSS, DES, Gaia and LSST data, we find that it is possible to measure the complete set of properties (including the phase-space coordinates during the flyby) of dark perturbers with , up to a 1d degeneracy between the mass and velocity.
keywords:cosmology: theory - dark matter - galaxies: haloes - galaxies: kinematics and dynamics - galaxies: structure - astrometry, surveys, proper motions
Embarrassingly few observational techniques exist with which to reliably gauge the presence of low-mass Dark Matter (DM) clumps in the Universe around us. Even though swarms of completely dark subhaloes are predicted to orbit galaxies like our own (Diemand et al., 2008; Springel et al., 2008), there is still no proven tactic to detect them. Until recently, the only acclaimed direct probe has been strong gravitational lensing (see e.g. Dalal & Kochanek, 2002; Hezaveh et al., 2014). However, such measurements are statistical in nature: only the total mass fraction locked in small sub-structure can be typically estimated. The largest of the small subhaloes, those with masses between and can perhaps be identified individually (see e.g. Vegetti et al., 2012; Hezaveh et al., 2013). Yet, there is little hope to prove that they do not possess a substantial baryonic component.
In the Milky Way, arguably, the frontier of low-luminosity low-mass substructure detection was breached long ago. A population of DM-dominated dwarf satellites is known to surround the Galaxy (see e.g. Aaronson, 1983; Mateo, 1998; Tolstoy, Hill & Tosi, 2009). For every solar mass of hydrogen, these dwarfs contain hundreds of solar masses of Dark Matter (see e.g. Walker et al., 2009). Granted, the total masses of the subhaloes the dwarf spheroidals inhibit are not currently known. The initial mass of the faintest of them can be as high as (see e.g. Wheeler et al., 2015). However, it is not implausible that they could have started with as little as of DM (Bland-Hawthorn, Sutherland & Webster, 2015). By redshift zero, their dark masses have probably been depleted further while their luminosities are at the level barely seeable at distances far smaller than cosmological (see e.g. Belokurov, 2013).
Between and , there ought to be roughly 2000 DM subhaloes in a host if DM were cold (Springel et al., 2008). As of today, the census of the Milky Way’s dwarf population is incomplete, but the total is predicted not to exceed several hundred satellites (see Tollerud et al., 2008; Koposov et al., 2008, 2009). It is clear that more than half of the substructure in this mass regime must remain invisible. Recently, Warm Dark Matter has been invoked as a silver bullet to make the problems associated with low-mass substructure go away (see e.g. Lovell et al., 2012, 2014). However, in line with the up-to-date Ly measurements, if the DM were warm, its mass ought to be greater than 3.3 keV at the 2 level (see e.g. Viel et al., 2013). This, in turn, leaves very little wriggle room to modify the substructure mass function above km s (see Schneider et al., 2014).
The next goal, therefore, is to discover subhaloes totally devoid of stars. The most promising avenue is the detection and characterization of the perturbations DM clumps induce on stellar streams in the Galaxy (e.g. Ibata et al., 2002; Johnston, Spergel & Haydn, 2002). Through a series of numerical experiments, it has now been established that dark subhaloes should interact with stellar streams sufficiently frequently, leaving small but observable signatures in their wake (e.g. Siegal-Gaskins & Valluri, 2008; Carlberg, 2009, 2012; Yoon, Johnston & Hogg, 2011; Ngan & Carlberg, 2014; Ngan et al., 2015; Carlberg, 2015). For example, during the lifetime of a Pal 5-like stream roughly 5 direct impacts with subhaloes with are expected and 20 impacts with subhaloes with (Yoon, Johnston & Hogg, 2011). There are several tens of tidal streams already discovered in the Milky Way (Newberg & Carlin, 2016) - plenty of potential carriers of signs of past interactions.
Following the treatment introduced in Carlberg (2013), we laid down the basic equations governing the formation of stream perturbations in Erkal & Belokurov (2015, Paper I). Paper I presents a simple picture of the three phases of stream gap evolution, and, most importantly, demonstrates that a measurement of the gap density profile alone is not sufficient to break the strong degeneracy between the perturber’s mass and size, the impact geometry and kinematics, and the time since flyby. However, if complementary phase-space measurements around the gap were collected, would the evidence assembled be ample to reconstruct what exactly happened during the run-in? This is the question we strive to answer in this Paper.
In this work we will simulate cold streams similar to that of Pal 5 (see e.g. Odenkirchen et al., 2003). We will consider the 6d phase-space structure of the gap and show that using this expanded set of observables, it is possible to infer the properties of the subhalo’s flyby, up to a 1d degeneracy. This is a degeneracy between the subhalo mass and the velocity of its flyby which is unavoidable since it is actually a degeneracy in the velocity kick imparted on the stream particles during the flyby. We will show that with proper motions and distances from Gaia, Dark Energy Survey (DES)-like photometry, and moderate resolution spectroscopy for radial velocities, it will be possible to detect and characterize dark subhaloes down to within 10 kpc. We will also discuss the mass range which can be probed by upcoming surveys like the Large Synoptic Survey Telescope (LSST) and give general expressions for the mass sensitivity as a function of the errors in each of the 6d phase-space coordinates.
Finally, we deliberate on the constraints that can be derived if only a subset of the observations can be secured. For example, for sufficiently distant streams it may not be possible to get proper motions with Gaia. We show that in general one only needs to measure the following: the shape of the gap on the sky, the density profile of the gap, and the line-of-sight velocity along the gap. Fortunately, these are the simplest measurements to make for targets across a large range of heliocentric distances, so this does not limit our ability to infer subhalo properties.
This paper is organized as follows. We give the setup and motivation in Section 2. In Section 3, we present analytic expressions for the 6d structure of the gap. Next, in Section 4, we use these results to infer the properties of a subhalo from a gap in an N-body simulation without accounting for observational errors. In Section 5, we discuss the information content of each observable and find the minimal set of observables which are needed to constrain subhalo properties. In Section 6, we infer the properties of subhaloes from gaps in simulations using observational errors of current and future surveys. Continuing with realistic errors, we present analytic expressions for the mass range accessible as a function of survey errors in Section 7. In Section 8 we discuss complications which will arise when extending this model to more realistic streams. Finally, we conclude in Section 9.
2 Setup and Motivation
The aim of this work is to show that it is possible to start with a gap in a tidal stream and infer the properties of the subhalo which created it. In order to study this problem, we will consider a simple setup with a stream on a circular orbit which is perturbed by the flyby of a Plummer sphere. This is identical to the setup of Paper I. In that work, we were primarily interested in the density structure of the gap and we showed that a measurement of the stream density profile near the gap does not uniquely constrain the properties of the subhalo which created the gap. In fact, it only provides constraints on three out of the seven parameters which describe the subhalo flyby. Thus, more constraints are needed. The results of Paper I also contain predictions for the entire 6d structure of the gap which we will exploit in this work. With this extra information, we will show that we can constrain the subhalo properties down to a 1d degeneracy between mass and velocity.
One of the main benefits of this simplistic model for gaps is that all of the predictions are analytic. With these expressions for the stream shape, we can predict the subsequent evolution of the stream after a subhalo flyby. For a given gap, the model can then be used to infer the properties of the subhalo which created it. Furthermore, the analytic expressions will allow us to probe the constraining power of each observable and we will discuss what degeneracies remain if only certain observations are possible. Finally, we will use these expressions to give constraints on what mass range is gaugeable given the errors of current and future surveys. Without these analytic predictions, a lot of this would not be possible and we would need to build intuition entirely through individual numerical examples.
Figure 1 displays the geometry of the subhalo-stream encounter as used in this work. The stream is on a circular orbit and we choose our coordinates to be centered on the point of closest approach along the stream with denoting the radial direction, denoting the direction along the stream, and denoting the direction perpendicular to the stream plane. This is identical to the setup in Paper I. We assume that the stream is orbiting in a spherical potential which is otherwise unconstrained. The stream is at a radius with a circular velocity . This model allows for any impact geometry. We denote the subhalo velocity as and the stream velocity as . At the point of closest approach, the vector between the stream and the subhalo must be perpendicular to the relative velocity . This relative velocity, along with the impact parameter , defines the subhalo’s trajectory up to an overall sign since the subhalo could pass on either side of the stream. We will encode this sign in the impact parameter, . Next we have the mass of the subhalo, , and its scale radius, . Finally, we have the time since impact, . For a given gap, we take the impact point to be the current position of the gap center integrated backwards in time by . Thus we have seven model parameters: three velocities, , the impact parameter, , the mass , the scale radius, , and the time since impact, . Following the convention of Paper I, we break the relative velocity between the stream and subhalo into a velocity along the stream, , and a perpendicular velocity of . The magnitude of the relative velocity is denoted by .
During the flyby, we assume that the stream is moving along a straight line. This is needed to compute the velocity kicks along the stream. Since the velocity kicks take place over a relatively small distance, several kpc at most for the most massive subhaloes we will consider, this is a reasonable assumption.
3 New Observables
The first observable we will consider is the density profile along the stream as in Paper I. In addition, we will now use the analytic model produced in Paper I to make predictions for the 6d phase space structure of the stream after the flyby, giving five new observables. First we have the spatial shape of the stream near the gap which gives two observables: 1) the distance to the stream as a function of the angle on the sky, vs (where is the angle along the stream) and 2) the shape of the stream on the sky, vs . Next we have the velocity structure of the gap which gives three observables: 3) the radial velocity as a function of the angle on the sky, vs , 4) the proper motion along the stream as a function of the angle on the sky, vs , and 5) the velocity perpendicular to the orbital plane as a function of the angle on the sky, vs . We note that these observables are what an observer in the center of the galaxy would measure. They do not account for the solar position or velocity. We will give analytic predictions of these observables below and then use them to infer the subhalo from several N-body simulations. Note that these predictions can only be expressed parametrically, i.e. to describe the shape of the gap on the sky ( vs ) at time , we express both and in terms of the angle from the point of closest approach, .
Crucially, the structure of the gap is controlled by the velocity kicks imparted on the stream during the flyby. The kicks in the , , and direction (radial, along the stream, and perpendicular to the orbital plane), as a function of distance from the point of closest approach were computed analytically for a Plummer sphere perturber in Paper I:
These formulae are slightly different than the ones in Paper I since we have re-expressed the kicks in terms of the three components of the flyby velocity instead of two velocity components and an angle. We have also re-expressed the distance from the point of closest approach, , in terms of the angle along the stream, i.e. . Next, we list the prediction for each of the 6d phase space coordinates.
First we have the angle along the stream. This is given in equation 16 of Paper I:
where denotes the initial distance between a point along the stream and the point of closest approach,
and the timescale is given by
Note that these expressions are slightly different from those in Paper I since we have re-written the angle in terms of ratios of , , and . The density is related to this angle since .
The radial distance of each particle in the stream is computed in equation 8 of Paper I:
where the velocity kicks, , are specified in (1).
The kicks in the direction tilt the orbital plane of the resulting orbits by an angle of in the direction, therefore the offset in stream latitude of the particles is given by
where which denotes the difference in stream latitude from the unperturbed stream track.
The radial velocity can be derived from a time derivative of (8) which gives
The velocity along the stream is given by
where is the angular velocity along the stream and is given in equation 12 of Paper I and is given in (8). Plugging in these expressions and taking only the leading order terms we get
We note that we only consider the leading order terms in and in this equation since the analysis in Paper I was done at leading order.
The velocity in the direction comes from a time derivative of the offset in the direction, i.e. the time derivative of (9):
By combining these six results, (2) and (8)-(13), we have a parametric model for the phase space structure of the stream at time in terms of the initial position along the stream, . In Figure 2-Figure 4, we show a comparison of the model presented in this section against the results of N-body simulations which we will describe in the next section. From these Figures, it is evident that the model can successfully reproduce the structure of the gap.
4 Inference without errors
Now that we have introduced the new observables, we will use a suite of N-body simulations to demonstrate that the model works. We will then run a set of MCMCs to show that the model can be used to constrain the properties of the subhalo which created the gap.
To generate the gaps we carry out three N-body simulations with the pure N-body part of gadget-3 which is closely related to gadget-2 (Springel, 2005).
The setup is similar to the realistic simulations in Paper I. We first generate a stream by taking a globular cluster modelled as a King profile with a mass of , a core radius of 13 pc, and . This progenitor is placed on a circular orbit with a radius of 10 kpc around an NFW potential (Navarro, Frenk & White, 1997) with , and kpc. The circular velocity at 10 kpc is km/s. The King profile is represented with particles which have a smoothing length of 1 pc. The cluster is evolved for 5 Gyr and in this time a cold stream with a total length of is produced. We take this stream and place a subhalo on an orbit to impact the center of the leading arm of the stream, which corresponds to a radius of 9.80 kpc, at 150 km/s perpendicular to the orbital plane. We then follow the evolution of the stream for an additional 5 Gyr. We have chosen a simple impact geometry in this example but the analytic model works for any impact geometry.
In our experiments, we use three different subhaloes to showcase the range of masses and scale radii for which we will be able to infer the subhalo properties. These are chosen to roughly lie in the center of the scale radius-mass relation in the public data release of the Via Lactea II simulation (Diemand et al., 2008). First, we have the largest subhalo with with pc. Next, we have an intermediate subhalo with with pc. Finally, we have the smallest subhalo with and pc. We note that while these subhaloes were chosen in the middle of the scale radius-mass relation from Via Lactea II, we have modelled them with Plummer spheres instead of NFW profiles. As a result, our Plummer spheres reach a maximum circular velocity at a smaller radius than if they were modelled with NFW profiles. If we instead compare the radius at which the maximum circular velocity of our subhaloes is reached to those in Via Lactea II, we find that our subhaloes are still within the scatter but are now closer to the edge of the relation rather than the center.
As we saw in Section 3, the amplitude of the oscillations induced by the subhalo’s passage is governed by the maximum velocity kick induced by the subhalo. For a direct impact, this maximum is proportional to which follows from (1). The properties of these subhaloes are arranged so that the size of the perturbation decreases by a factor of two as we go from the most massive to the intermediate subhalo, and another factor of two as we go from the intermediate to the least massive subhalo. Finally, we have to choose how long after impact to look at each of the gaps. We choose to look at the perturbations when the deflection on the sky ( vs ) is relatively large at Myr. Note that the different subhaloes will have created different sized gaps by this time. While the gap will grow in time, becoming easier to observe, we choose to look at the gap after a short time where our model is most accurate. At late times, the distribution of energy and angular momenta along the stream can cause the gap to fill in which is not captured by our model. We choose to look at all subhaloes at the same time to make it a fair comparison.
In Figure 2,3,4 we show the six observables for each of the three instances of the perturbed stream at these times. Note that the scales on the axes are different in each of the figures. The grey bands in these figures show the width of the stream and the black curves show the prediction from the model for each setup. We stress that the black curves are not fits to the data but rather a prediction of the model using the true parameters. While the model agrees well with the simulations, it does not capture the entire behavior. For example, if we look at the panel in Figure 2 or Figure 3, we see that the model does not pass through the center of the grey band. This slight disagreement is important as it will result in biases when we fit our model to the simulations as demonstrated in Section 4 and Section 6.
In these simulations, we have focused on direct impacts which will produce the largest signatures. The analytic results of Section 3 allow us to explore how large the signatures will be if we change the impact parameter, the mass of the subhalo, or the velocity of the flyby. The size of the perturbation is governed by the size of the velocity kick which can also be estimated qualitatively by taking the typical acceleration during the flyby, , multiplied by the duration of the flyby, , to give a velocity kick of . As we showed in Figure 3 of Paper I, the maximum kick in the direction is given by so the qualitative approach almost gives the correct behavior. Thus we see that as we increase the mass of the subhalo, we will get a larger signature. Likewise, if we increase the velocity of the flyby, the impact parameter, or the subhalo’s size, we will get a smaller signature. The maximum size of the velocity kicks in the and can also be computed using (1) and exhibit roughly the same behavior although the dependence on the velocities, , and is more complicated.
We will now infer the properties of the subhalo without accounting for any observational error. This showcases how well this method works in the best case scenario. In Section 6 we will do the same exercise accounting for observational errors for both contemporary and future surveys. We use emcee (Foreman-Mackey et al., 2013) to explore the likelihood landscape of our models.
The MCMC samples the seven model parameters, , as well as 6 nuisance parameters for overall shifts in the 6 phase space coordinates. These nuisance parameters are needed to account for possible offsets of the center of the gap due to epicyclic motion. As we can see from (1), the center of the gap receives no kick along the stream so it has the same orbital period as the unperturbed stream. However, it can receive kicks in the other directions which causes the center of the gap to execute epicycles around its unperturbed position. As a result, when we measure the location of the gap, we do not know the current offset of the center of the gap from the unperturbed stream and we need to marginalize over this.
The model described in Section 3 gives the centroid of the gap in all six observables. The centroid of the gap in the N-body simulations is found by fitting a Gaussian profile (controlled by 2 parameters describing its center and intrinsic width) to each observable in 1 degree bins along the stream. This gives the centroid as well as its error for all observables except the density which is given Poisson errors. This procedure mimics the process of decoupling the stream and the background, where, at least for some of the observables, the separation is done in bins of angle along the stream and not per star.
With regards to modelling the subhalo-stream encounter itself, the likelihood for each observable is simply:
where denotes the angle bins along the stream, are the centroid and its error in the N-body simulation, and is the centroid from the model. We assume the observables are independent so the full log-likelihood is the sum of the individual log-likelihoods for each observable.
In the analysis presented here, our priors are as follows. For the mass, we assume a uniform prior in between a mass of and . For the scale radius, we assume a uniform prior between pc and kpc. For the impact parameter, we take a uniform prior between kpc and kpc. The prior on the velocity of the subhalo is a Maxwellian distribution with a dispersion, , related to the circular velocity at kpc, , by the relation . For our setup this gives km/s. This prior is very important since, as we will see and discuss below, there is a degeneracy between the mass and velocity of the flyby. For the time since impact, we assume a uniform prior between Myr and Gyr. Next we have priors for the nuisance parameters. The shift in the angle along the stream is uniform from to . The shift in the radial distance is uniform from kpc to kpc. The shift in all velocity components is uniform from to km/s.
In Figure 5 we show the result of the MCMC modelling for the lowest mass subhalo. We see that without accounting for observational errors, this method can accurately infer the properties of subhalo with a mass of . Almost all of the parameters are well constrained. The only exception is the degeneracy between the mass and the velocity which we will discuss in detail in Section 5. This degeneracy is due to the fact that the model cannot infer and independently, but rather can only deduce their ratio . As a result, the constraint on the mass depends on the prior we place on the subhalo velocity. Note that the likelihood does not peak at the true value for every parameter. In fact, and especially are discrepant. This is likely due to the fact that our model does not account for all of the features in the stream, such as epicyclic feathering (e.g. Küpper, MacLeod & Heggie, 2008; Küpper et al., 2010) and that we assume that the stream is on a circular orbit. While the progenitor itself is on a circular orbit, the stream it generates will be slightly eccentric. For the subhaloes with and , the figures look similar albeit with slightly less scatter.
5 Degeneracy and Information Content
As we saw in Figure 5, there is a degeneracy between the mass and velocity of the subhalo. We will now study the set of degeneracies that arise when inferring stream properties.
5.1 Degeneracy in the velocity kicks
First, we will study the degeneracy present in the velocity kicks, , induced by the subhalo’s passage. These kicks govern the subsequent evolution of the gap so the best we can do is recover the information content encoded in the kicks. These kicks are given in (1), which we will re-write as
where we have re-expressed (1) in terms of the six coefficients, , which depend only on the parameters , which describe the flyby. For completeness we give the in (42). The entire effect of the flyby is encapsulated in these three kicks which are a function of . If any of these coefficients are different, a different gap will be produced. In the ideal case, we would measure these functions directly and then read off the coefficients. Degeneracies can only arise if there are directions we can move in the 6d parameter space, , which do not change any of these coefficients. If we consider a small perturbation in our parameter space, , the degeneracies would need to satisfy a system of linear equations:
where the gradient is taken with respect to . This is a system of 6 equations with 6 unknowns. In principle, this could have one solution and no degeneracy. However, for generic parameters, these 6 constraints are not linearly independent and give rise to only 5 independent constraints. Therefore, there is a 1d degeneracy which involves a scaling of mass and velocity:
This degeneracy is insurmountable since it does not affect the kicks themselves and hence has no effect on the gap.
5.2 Degeneracies in the observables
We can now apply this same argument to the observables described in Section 3. We can then combine various subsets of these observations to see what these combinations can constrain. The best we can hope for is to match the ideal case of a 1d degeneracy. The observables have an additional parameter which must be fit: the time since impact, . Thus when considering observables we have a 7d parameter space which we will denote, .
Next we have the radial distance along the stream, (8), which we re-write as
The stream latitude is governed by (9), which we re-write as
The radial velocity is given by (10), which we re-write as
The proper motion along the stream, , is given in (12), which we re-write as
The proper motion in the direction is give in (13), which we re-write as
Note that the denominator is identical for all of these observables. This is because all of these terms are related to the velocity kicks imparted by the subhalo which have the same denominator (see the denominator of (1)). For completeness we give the in (43).
In the expressions for the observables above, (18) - (23), there are 13 separate coefficients which govern the observables. If all observables are measured, we can determine the resulting degeneracies in the parameters, , by trying to find variations of the parameters, i.e. , which do not change the :
where the gradient is taken with respect to . By carrying out this procedure, we find that these 13 constraints on 7 unknowns give rise to 6 linearly independent constraints, leaving the same degeneracy as above, i.e. (17). This means that if all observables are measured, we can recover the same information as if we measured the kicks themselves.
Note that this result is for a generic value of the parameters. For specific values of the parameters, it is possible to have additional degeneracies if some of the coefficients vanish, i.e. if for some . For example, since the observables oscillate in time, there will be times when some of the observables are flat and contain less information. Another example is if . In this case, as we can see from the velocity kicks in (1), and only depend on the combination and not or independently. As a result, if the observables do not constrain (or poorly constrain) the shape in and , we should expect an additional degeneracy to appear between and where is constant. This is not a true degeneracy since the gap size does not only depend on . However, since the leading order term in the gap size only depends on , the gap size is not very sensitive to . As a result, and can be effectively degenerate. We will see an example of this in Section 6.4.
5.3 Degeneracies in typical observational setups
Now that we have discussed how to classify the degeneracy when using all observables, we will consider the degeneracies present in common observational setups which only use certain observables. We will repeat the procedure above, using only the coefficients present in the selected observables.
5.3.1 Gap density, gap shape in , and radial velocity
One of the simplest useful setups involves measuring the density of the stream on the sky, the shape of the stream on the sky, and the radial velocity profile along the stream. These three observables depend on the coefficients , and . Using these coefficients, we assemble and solve the system of 7 linear equations as in (24). Interestingly, these 7 equations also result in 6 linearly independent constraints on the 7d parameter space, , leaving the degeneracy between mass and velocity, i.e. (17). This setup is quite simple since it only requires the radial component of the velocity which is far easier to measure than proper motion. This approach should also allow us to extend the search for gaps out to larger distances since the errors for these observables do not depend as sensitively on distance as the errors in proper motion.
5.3.2 Gap density and gap shape in
Another approach is to measure the gap density and gap shape on the sky which can be computed using only photometric surveys and do not require any velocity information. These observables depend on the coefficients , and . We setup and solve the 5 linear equations as in (24) and find 5 independent constraints, leaving a 2d degeneracy. One of the dimensions of this degeneracy is the mass-velocity degeneracy discussed above. The second degeneracy is more complicated and is between all of the parameters.
5.3.3 Gap shape, radial velocity, proper motions
In this approach, we make use of all of the observables except the density along the stream, i.e. , , , , and vs . This setup makes use of all the coefficients except and . We setup and solve the 11 linear equations as in (24) and find 6 independent constraints leaving the degeneracy between mass and velocity. This is a useful setup since we need to understand the stripping rate of the progenitor as a function of time to fully model the density along the stream. Since this may prove challenging, it is useful to know that the density along the stream is not needed to constrain the subhalo properties.
In Section 6.4 we will compare the inferred subhalo properties when using these three subsets of observables against the inferrence when using all observables.
6 Inference with imperfect data
Now we will repeat the inference done in Section 4.2 and account for observational errors.
6.1 Observational Errors for Current and Future Surveys
We will now assume three different observational setups to estimate the subhalo mass range accessible with this technique. Their details are available in Table 1. The first setup we consider follows closely the properties of the SDSS (Ahn et al., 2012) data. Next, we have the “Super-SDSS” survey which consists of DES-like (The Dark Energy Survey Collaboration, 2005) photometry, coupled with Gaia (Perryman et al., 2001) astrometry for proper motions, and deep multi-object spectroscopy from VLT (e.g. Koposov et al., 2011; Gilmore et al., 2012), WEAVE (Dalton et al., 2012), and 4MOST (de Jong et al., 2012) for radial velocities. Finally, we have the LSST setup, which is similar to Super-SDSS but uses LSST (LSST Science Collaboration et al., 2009) photometry and Gaia+LSST for proper motions to extend observations to much fainter stars.
|SDSS||0.1”,||5 km/s at , 20km/s at||None||None|
|Super-SDSS||0.1”,||1 km/s at , 5 km/s at||Gaia||DES photometric accuracy|
|LSST||0.1”,||1 km/s at , 5 km/s at||Gaia+LSST||LSST photometric accuracy|
6.2 Realistic Stellar Observables
In order to understand the subhalo mass range which can be probed by these surveys, we need to populate the stream with a realistic distribution of stars. We use the N-body simulations described in Section 4.1 which model the globular cluster with identical mass particles. We use PARSEC stellar evolution models (Bressan et al., 2012; Chen et al., 2014; Tang et al., 2014) to assign masses and luminosities to each particle. We chose the initial mass function (IMF) from Chabrier (2001) and used a single isochrone with an age of 10 Gyr, and a metallicity of 0.01 . We made use of two photometric systems: SDSS and Gaia . The stellar mass of each particle was sampled from the corresponding present day mass function until the total stellar mass was equal to , the mass of the cluster used in the simulation. This resulted in approximately 360,000 N-body particles being tagged as stars. For each star, the magnitude in each band was then linearly interpolated from the tabulated values in the isochrone table. Given the magnitudes of each particle, observational noise was added to the position, velocity, and magnitudes using the errors and cuts in Table 1. We note that the number of stars above the -band cutoff in Table 1 and within 20 degrees of the gap center increases dramatically as we consider the deeper surveys. The SDSS setup has 1404 stars, the Super-SDSS setup has 4039 stars, and the LSST setup has 10735 stars. We stress that we do not account for contamination from Milky Way stars or from background galaxies.
The relative distances along the stream were calculated using the photometric parallax method. Given the isochrone of the stream, we use the color to infer the r-band absolute magnitude. The apparent -band magnitude combined with the absolute magnitude gives the distance modulus and hence the distance. We propagate the errors from the photometric accuracy for each color and use the slope of the isochrone to give the error in absolute magnitude. We do not account for uncertainties in the isochrone itself as we are primarily interested in relative distance changes across the gap. For the SDSS setup, we chose not to use any distance measurements since the distance errors are too large to give useful constraints. For the Super-SDSS setup, we use the expected photometric accuracy of DES. To estimate the single epoch photometric accuracy, we use the photometric catalog described in Koposov et al. (2015) which is derived from publicly available images from the first year of DES. We assume that all five epochs of DES will have a similar accuracy to get a final photometric accuracy which is smaller than the single epoch accuracy. For completeness, we list our expected five epoch DES -band magnitude errors: 0.0005, 0.0096, 0.0019, 0.0038, 0.0079, 0.017, 0.039, 0.091 for -band magnitudes of 17,18,19,20,21,22,23,24 respectively. The -band magnitude errors are similar to these. For the LSST setup, we use the LSST photometric accuracies which are reported in LSST Science Collaboration et al. (2009).
6.3 MCMC Details
The details of the MCMC are almost identical to those presented in Section 4.2. The only difference is that the observables for each star are now scattered by the appropriate observational error and this is accounted for when we compute the centroid and its error for each observable along the stream.
In Figure 6,7,8 we show the six observables and their error bars for the subhalo for the three different observational setups: LSST, Super-SDSS, and SDSS. We note that the error bars in these figures are the errors on the mean within each angular bin and are therefore significantly smaller than the error on any individual measurement. We see that, as expected, the shape of the stream on the sky and the density have the smallest error bars, followed by the distance along the stream, and then the radial velocity. We see that the error bars for the Super-SDSS setup are already quite constraining. As we move from Super-SDSS to LSST, the main improvement is the superior photometry and depth (i.e. number of stars) which in combination give smaller error bars in the shape of the stream gap on the sky as well as its density. There is also a noticeable improvement in the proper motion measurements.
Figure 9 shows the 1 constraints on and for each subhalo and each observational setup. For and we see that the constraints get tighter as we use higher quality data. The exception is where the constraints do not improve with better data. This is due to the degeneracy between and which means that the constraint on comes mostly from the prior on the subhalo velocity.
Figure 10 gives the 1 constraints on and for the subhalo using LSST errors. In this figure, we vary the observables used in the MCMC. We see that if we use all of the information, or subsets of the observables which can constrain all of the properties, we obtain results of similar quality. However, if we have too few constraints then the error bars get larger. The three subsets we use here are described in Section 5.3.
In Table 2 we give the 5-95% confidence intervals for and for each subhalo and observational setup. The trends here are identical to those seen in Figure 9. In Table 3 we show the same intervals for the subhalo when we only use , , and density information. As we discussed in Section 5, these observables are sufficient to constrain the subhalo properties up to the 1d degeneracy between mass and velocity. As expected, the constraints are similar to when we use all of the information.
|Survey||, pc||, pc||, pc|
Figure 11 presents the posterior distribution for the subhalo when the LSST errors are used. This figure can be compared against Figure 5 where we showed the posteriors with no observational error. While there is more scatter for the LSST errors, the overall trend is similar. Most of the parameters are well constrained except for the degeneracy between and . There is a hint of the degeneracy between and where is constant and which we mentioned in Section 5. This can best be seen in the 2d posteriors of and . This degeneracy will be much more vivid when we use SDSS errors.
Figure 12 demonstrates the posterior distribution for the subhalo when the SDSS errors are used. This setup has significantly more scatter than the LSST setup. The degeneracy between and is present as expected. In addition, the degeneracy between and when is now evident. It is especially clear in the 2d posterior of . This degeneracy is due to the fact that the SDSS setup has no leverage on the distance to the stream and only a modest constraint on the radial velocity. As we discussed in Section 5, this allows for an additional degeneracy to appear. Note that this is not a true degeneracy since we have some constraint on the radial velocity and the gap size has a weak dependence on . As a result, the likelihood is still peaked near the true value of .
7 Mass Range Probed by Current and Future Surveys
In Section 6, we showed the results of a suite of MCMCs where we were able to measure the properties of subhaloes down to . While this explicit demonstration is useful, we can also use the analytic results of Section 3 to build some intuition as to the constraining power of each observable. We will restrict the analysis to the case of a direct impact, i.e. , since these create the largest signature. The detectability of a gap is controlled by the maximum kick imparted by the subhalo. This maximum kick is given by the maxima of (1):
Thus we see that the maximum velocity kicks are roughly
We can now plug this expression into those for each of our observables and compare the maximum deviations of the stream against a typical uncertainty on the centroid from the LSST setup. We note that these errors are significantly less than the measurement error for a single star.
The maximum radial deflection is given by the maximum of (8) which gives
This must be larger than the error in estimating the stream centroid, . To determine the mass range accessible with a given error, we plug in the expression for , i.e. (26), and get
where is the relative velocity of the flyby so we see that lower velocity flybys allow us to probe a lower mass range. We have re-arranged this expression to isolate the parameters describing the subhalo since these are the most relevant for understanding Galactic substructure. For typical values of km/s, kpc, (i.e. a flat rotation curve), and a radial distance error of pc, we get
This expression can be used to understand the mass range probed as we change the physical and observational setup. At the edge of this limit, i.e. if , we should expect a signal-to-noise ratio of 1. Note that we have not included the dependence in this expression.
The maximum deflection perpendicular to the plane is given by
This must be larger than the error in estimating the stream’s position on the sky, . This gives
For the typical values of km/s, with the error in the stream latitude is given by the error when estimating the centroid, , we get
The maximum deflection in radial velocity is given by
This must be larger than , which gives
For typical values with km/s, , and km/s, this gives
The maximum deflection in proper motion along the streams’s orbit is given by
This must be larger than , giving
For typical values with km/s, , and km/s at 10 kpc, this gives
The maximum deflection in the proper motion perpendicular to the orbital plane is given by
This must be larger than , which gives
For typical values of km/s, , and km/s at 10 kpc, this gives
By comparing these results, we see that the most discerning observable is the shape of the stream on the sky, followed by the distance to the stream, the radial velocity, and finally the proper motions. The proper motion perpendicular to the orbital plane is less constraining since the velocity oscillations are larger along the orbit. It is not possible to make such a simple estimate for the constraining power of the density along the stream. Note that these formulae can be used to quickly evaluate the mass and size range probed by any survey.
The results above show that with current and upcoming surveys, it is possible to infer the properties of subhaloes down to a mass of at a distance of 10 kpc (note, however, that according to Section 5.3, this measurement can be extended to much larger distances). In addition, we quantify the accuracy with which the geometry and the dynamics of the encounter as well as the mass and the size of the perturber can be measured. Importantly, we demonstrate that while there exists a degeneracy between subhalo’s mass and velocity, its detrimental effect can be remedied with the use of reasonable priors. We also presented expressions for the mass range accessible for a given observational error in Section 7.
These mass limits should be considered as a best case scenario. In reality, there will be several additional complications which we will now discuss. These will introduce additional errors which may raise the mass limit measurable with the method discussed here. We will also discuss other causes of substructure in tidal streams and how these can be distinguished from gaps caused by subhaloes.
8.1 Discrete Degeneracies
In Section 5, we characterized the degeneracies when inferring subhalo properties. We showed that if we use all of the observables, there is a continuous degeneracy between and where is fixed. In addition to this continuous degeneracy, it is possible that there are discrete degeneracies with different points in the parameter space giving rise to the same gap shape. Since the gap shape oscillates in time, if we wait roughly an orbital period, the gap shape will look similar in vs . Of course, the gap will now be larger so we must change the other parameters to get the same gap size and shape. Thus, we might imagine a discrete degeneracy is possible.
We have searched for this discrete degeneracy by running gradient descents in our negative likelihood space and initiating MCMCs near the location of different likelihood peaks. In the case of no observational errors, we have found no evidence for such degeneracy since the log-likelihoods of the varying peaks differ by more than 1000 from the true peak. As a result, we do not believe this is a true degeneracy. However, if we include observational errors, discrete degeneracies are possible in some cases. For the subhalo and using SDSS observations, we cannot distinguish between solutions which differ by an orbital period. However, if we use Super-SDSS or LSST observations, these solutions differ by roughly 5 and 10 respectively in log-likelihood. For the less massive subhaloes, i.e. and , the discrete degeneracies cannot be distinguished even with LSST quality observations. In this work we have only considered one setup: a direct impact on the stream which we observe when the oscillation in the shape on the sky is large. It is not clear how these effective degeneracies will change for a different impact geometry, if we observed the gap at a different time, or if the progenitor is on an eccentric orbit. We will explore this issue in future work.
8.2 Eccentric orbits
The analysis done in this work was for progenitors on circular orbits. In reality, progenitors will be on orbits with a range of eccentricities. This will result in the gap stretching and compressing as it passes near pericenter and apocenter respectively. Since analytic orbits are not known for an arbitrary spherical potential, it is not possible to repeat the analysis presented here for eccentric orbits. However, the procedure outlined in this work can be done numerically and the qualitative picture presented here should be the same. In future work, we plan to address this by re-casting this analysis in action-angle formalism.
8.3 Uncertainty in the potential
In the method above, we have assumed perfect knowledge of the potential. Any uncertainty in the potential, and thus the shape of the stream, will complicate the inference of the subhalo properties. However, since the potential is expected to be smooth, its effect on the uncertainty of properties derived from kinks and wiggles over scales in the stream should not be very large. These wiggles can only be created by small scale structure.
One additional complication which we have mentioned above is the fact that all of these observables oscillate in time. Thus, the quality of the inference depends on what phase of the oscillation the gap is in. In the worst case scenario, the observables which have the most constraining powers (i.e. the shape on the sky or the radial velocity) will give less information than expected since they are in a phase with little signal. This will lead to a broader posterior and may change the results of Section 7 so that, in a particular phase, the proper motion observables are more informative than the shape of the gap on the sky.
8.5 NFW Subhaloes
In this analysis, we have assumed that subhaloes can be modelled as Plummer spheres. In reality, the subhaloes will be closer to NFW profiles (Navarro, Frenk & White, 1997). Since the kick induced by the passage of an NFW subhalo cannot be computed analytically, this means we will have to compute the kicks numerically. However, since numerics will already be necessary because of the non-circular orbits, this should not introduce much additional complication. Furthermore, since the NFW profile is controlled by two parameters (concentration and virial mass), this more realistic setup will have the same dimensionality as in this work and should exhibit the same qualitative features seen in this work.
8.6 Epicyclic Feathering
One physical process can produce density variations in a stream is epicyclic feathering (e.g. Küpper, MacLeod & Heggie, 2008; Küpper et al., 2010). In order to carry out the process described in this work on real streams, we will need to start with a detailed stream model. This model will include these epicyclic feathers so we will be able to determine the likelihood of a gap feature having arisen from the epicyclic feathering. More generally, the epicyclic feathers have different features from the gaps: the epicycles resemble cycloids while the stream near a gap looks like a sinusoid. Thus, with sufficient observations, it should be possible to distinguish the two. This will be easiest if the gap is far from the progenitor, probing an older part of the stream where it is well-ordered. Fortunately, these older parts of the stream are also the most likely to have gaps since they have been around the longest. On a promising note, the N-body simulations in our work naturally included this epicyclic feathering and our model was robust enough to make accurate predictions despite not accounting for these effects.
8.7 Gap structure from other substructure
In addition to dark substructure in the Galaxy, there are also baryonic structures like giant molecular clouds and spiral waves. These can potentially cause confusion but as we have shown in Section 4 and Section 6, the time since impact can be well constrained. This allows us to rewind the position of the gap back to where the impact occurred. If this results in the gap being created in or very close to the disk, baryonic effects should be considered. However, if the impact occurs far from the disk, this should allow us to rule out such effects. In addition, the results of this work show it is possible to get independent constraints on the mass and scale radius. This should give even more evidence to rule out baryonic structures if the gap is created from a large and low-density structure.
8.8 Searching for Stream Gaps
The primary signatures of gaps are the oscillations in the stream on the sky outlined in this work and the contrast between the low central density and enhanced density on the edge of the gap discussed in Paper I. The results of this work and the results of Paper I can be used to create match filters to search for these features. However, a simpler approach may be to search for qualitative features in the data. Namely, an underdensity with overdensities around it, and an oscillation in the shape of the stream on the sky. This requires only positions of stars from photometric survey data and does not require any velocities. Once such a feature is found, follow-up would be needed to constrain the radial velocities. As we discussed in Section 5.3.1, these three measurements would be sufficient to constrain the subhalo properties.
In this work, we have shown for the first time, how the properties of a dark subhalo can be extracted from the measurements of the gap it creates in a stellar tidal stream. We used the analytic results of Paper I to forecast the evolution of the 6d structure of the gap shape arising from a subhalo flyby. We then used these predictions to model gaps produced in N-body simulations. We carried out inference in an error-free case as well as for observables with realistic noise for a range of plausible instrumental setups. The main conclusions of our study are as follows:
(i) The information content of the subhalo-stream interaction is encoded in the velocity kicks received by the stream stars.
(ii) Because the velocity kick is an integral of the acceleration imparted by the subhalo onto stream particles, it depends both on the strength of the perturber’s pull as well as the duration of the encounter. In other words, there exists a degeneracy between the mass of the subhalo and the relative velocity of the flyby.
(iii) Using the shape of the gap without observational errors, it is possible to uniquely infer the impact geometry, its dynamics as well as the age for subhaloes with masses .
(iv) Importantly, it is also possible to independently infer the properties of the dark perturber, such as its scale length and its mass.
(v) If we account for observational errors from ongoing and planned surveys, it is possible to infer the impact geometry, its dynamics as well as the age for subhaloes with masses . For example, for the subhalo with a scale radius of 250 pc, the 1 uncertainty on the measurement of the scale radius, the impact parameter, the time since flyby decreases from 73pc, 135pc, 34.8 Myr for SDSS to 24 pc, 56 pc and 10.7 Myr for LSST.
(vi) Due to the mass-velocity degeneracy mentioned above, while it is possible to infer the mass of the perturber, the accuracy of the measurement does not scale simply with the signal-to-noise ratio of the observables. Yet, in our experiments we found the uncertainty of the perturber mass to be typically on the order of 0.15 dex.
(vii) Since many of the gap properties oscillate in time, a very similar gap will be produced by subhalo hits separated by an orbital period. While these different solutions can be ruled out in the error-free case, if we add observational errors these multiple solutions cannot be distinguished for the and subhaloes.
(viii) Using the analytic model for the stream, we explore the information contained in several different combinations of observables. We find that just three measurements, the stream density, the stream shape on the sky, and the radial velocity along the stream contain enough information to constrain the subhalo properties. We believe, therefore, that the stream gap analysis presented here can be carried out at heliocentric distances much larger than presented here.
The analytic nature of the results in this work allows us to build a simple intuitive picture of how to infer subhalo properties from stream gaps. We discussed the difficulties which could arise when extending this model to gaps created in more realistic streams. This will require numerical modelling but we expect the same qualitative picture we have described here. It should be possible to constrain subhaloes which impact streams on non-circular orbits up to the same 1d degeneracy seen here. In future work, we will explore this by re-casting this framework in angle-action coordinates.
Given the promising results of this study, the outlook for dark subhalo detection in the Milky Way is bright. If Dark Matter exists, we expect that these elusive clumps will soon be unveiled.
We would like to thank the Streams group at Cambridge for stimulating discussions and in particular we thank Sergey Koposov for help in understanding emcee. We thank Jason Sanders and Jo Bovy for useful discussions and we thank Wyn Evans for comments on an early version of the manuscript. We note that all of the MCMC triangle figures were made using triangle.py (Foreman-Mackey et al., 2014) The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 308024. We thank the referee for their thorough reading of our manuscript and their comments which improved the clarity of this work.
- Aaronson (1983) Aaronson M., 1983, ApJ, 266, L11
- Ahn et al. (2012) Ahn C. P. et al., 2012, ApJS, 203, 21
- Belokurov (2013) Belokurov V., 2013, NewAR, 57, 100
- Bland-Hawthorn, Sutherland & Webster (2015) Bland-Hawthorn J., Sutherland R., Webster D., 2015, ApJ, 807, 154
- Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
- Carlberg (2009) Carlberg R. G., 2009, ApJ, 705, L223
- Carlberg (2012) Carlberg R. G., 2012, ApJ, 748, 20
- Carlberg (2013) Carlberg R. G., 2013, ApJ, 775, 90
- Carlberg (2015) Carlberg R. G., 2015, ArXiv e-prints
- Chabrier (2001) Chabrier G., 2001, ApJ, 554, 1274
- Chen et al. (2014) Chen Y., Girardi L., Bressan A., Marigo P., Barbieri M., Kong X., 2014, MNRAS, 444, 2525
- Dalal & Kochanek (2002) Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
- Dalton et al. (2012) Dalton G. et al., 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 0
- de Jong et al. (2012) de Jong R. S. et al., 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 0
- Diemand et al. (2008) Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J., 2008, Nature, 454, 735
- Erkal & Belokurov (2015) Erkal D., Belokurov V., 2015, MNRAS, 450, 1136
- Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
- Foreman-Mackey et al. (2014) Foreman-Mackey D., Price-Whelan A., Ryan G., Emily, Smith M., Barbary K., Hogg D. W., Brewer B. J., 2014, triangle.py v0.1.1
- Gilmore et al. (2012) Gilmore G. et al., 2012, The Messenger, 147, 25
- Hezaveh et al. (2014) Hezaveh Y., Dalal N., Holder G., Kisner T., Kuhlen M., Perreault Levasseur L., 2014, ApJ, submitted; arXiv:1403.2720
- Hezaveh et al. (2013) Hezaveh Y., Dalal N., Holder G., Kuhlen M., Marrone D., Murray N., Vieira J., 2013, ApJ, 767, 9
- Ibata et al. (2002) Ibata R. A., Lewis G. F., Irwin M. J., Quinn T., 2002, MNRAS, 332, 915
- Johnston, Spergel & Haydn (2002) Johnston K. V., Spergel D. N., Haydn C., 2002, ApJ, 570, 656
- Koposov et al. (2008) Koposov S. et al., 2008, ApJ, 686, 279
- Koposov et al. (2015) Koposov S. E., Belokurov V., Torrealba G., Evans N. W., 2015, ApJ, 805, 130
- Koposov et al. (2011) Koposov S. E. et al., 2011, ApJ, 736, 146
- Koposov et al. (2009) Koposov S. E., Yoo J., Rix H.-W., Weinberg D. H., Macciò A. V., Escudé J. M., 2009, ApJ, 696, 2179
- Küpper et al. (2010) Küpper A. H. W., Kroupa P., Baumgardt H., Heggie D. C., 2010, MNRAS, 401, 105
- Küpper, MacLeod & Heggie (2008) Küpper A. H. W., MacLeod A., Heggie D. C., 2008, MNRAS, 387, 1248
- Lovell et al. (2012) Lovell M. R. et al., 2012, MNRAS, 420, 2318
- Lovell et al. (2014) Lovell M. R., Frenk C. S., Eke V. R., Jenkins A., Gao L., Theuns T., 2014, MNRAS, 439, 300
- LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv:0912.0201
- Mateo (1998) Mateo M. L., 1998, ARA&A, 36, 435
- Navarro, Frenk & White (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
- Newberg & Carlin (2016) Newberg H. J., Carlin J. L., eds., 2016, Astrophysics and Space Science Library, Vol. 420, Tidal Streams in the Local Group and Beyond: Observations and Implications. Springel
- Ngan et al. (2015) Ngan W., Bozek B., Carlberg R. G., Wyse R. F. G., Szalay A. S., Madau P., 2015, ApJ, 803, 75
- Ngan & Carlberg (2014) Ngan W. H. W., Carlberg R. G., 2014, ApJ, 788, 181
- Odenkirchen et al. (2003) Odenkirchen M. et al., 2003, AJ, 126, 2385
- Perryman et al. (2001) Perryman M. A. C. et al., 2001, A&A, 369, 339
- Schneider et al. (2014) Schneider A., Anderhalden D., Macciò A. V., Diemand J., 2014, MNRAS, 441, L6
- Siegal-Gaskins & Valluri (2008) Siegal-Gaskins J. M., Valluri M., 2008, ApJ, 681, 40
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Springel et al. (2008) Springel V. et al., 2008, MNRAS, 391, 1685
- Tang et al. (2014) Tang J., Bressan A., Rosenfield P., Slemer A., Marigo P., Girardi L., Bianchi L., 2014, MNRAS, 445, 4287
- The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration, 2005, arXiv:astro-ph/0510346
- Tollerud et al. (2008) Tollerud E. J., Bullock J. S., Strigari L. E., Willman B., 2008, ApJ, 688, 277
- Tolstoy, Hill & Tosi (2009) Tolstoy E., Hill V., Tosi M., 2009, ARA&A, 47, 371
- Vegetti et al. (2012) Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341
- Viel et al. (2013) Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, PhysRevD, 88, 043502
- Walker et al. (2009) Walker M. G., Mateo M., Olszewski E. W., Peñarrubia J., Wyn Evans N., Gilmore G., 2009, ApJ, 704, 1274
- Wheeler et al. (2015) Wheeler C., Oñorbe J., Bullock J. S., Boylan-Kolchin M., Elbert O. D., Garrison-Kimmel S., Hopkins P. F., Kereš D., 2015, MNRAS, 453, 1305
- Yoon, Johnston & Hogg (2011) Yoon J. H., Johnston K. V., Hogg D. W., 2011, ApJ, 731, 58
Appendix A Coefficients