© ESO 2018

1. Introduction

The determination of the orbital properties of the dwarf galaxies surrounding the Milky Way (MW) is a crucial step for unraveling the formation and evolutionary paths of these galaxies. Did the dwarf spheroidals (dSphs) and ultra-faint dwarfs (UFDs) lose their gas due to internal stellar feedback or due to a combination of internal and external effects, such as UV-background heating and/or tidal and ram-pressure stripping from the large Local Group spirals (e.g., Mayer et al. 2006; Bermejo-Climent et al. 2018; Revaz & Jablonka 2018)? Are the metallicity gradients observed for several of the Milky Way “classical” dwarf spheroidal galaxies (e.g., Tolstoy et al. 2004; Battaglia et al. 2006, 2011) an intrinsic property of these systems or could interactions with the Milky Way have had a role in setting their presence or erasing them in some cases (Sales et al. 2010)?

From a theoretical perspective, depending on the orbital and internal characteristics of the object (e.g., in terms of the mass density profile of its dark matter halo) it has been argued that such interactions could potentially transform initially disky/rotating dwarfs in spheroidal/pressure-supported systems (e.g., Mayer et al. 2001, 2006; Kazantzidis et al. 2017); however, the existence of very isolated dwarf galaxies, or those likely on their first approach towards the MW, having a spheroidal morphology and being in the process of losing/exhausting their gas would suggest otherwise (e.g., VV 124, Leo T, Phoenix I). Nonetheless, signs of tidal disturbance have been detected in the spatial distribution of the stellar component of some classical dwarf galaxies and ultra-faint systems (e.g., Carina I, Hercules I, Boötes I, see, e.g., Battaglia et al. 2012; Roderick et al. 2015, 2016) and it is still debated whether the presence of tidally unbound stars in spectroscopic samples can alter inferences of the dark matter (DM) distribution in these galaxies obtained from the observed l.o.s. velocity dispersion, in particular in the outer parts (e.g., ᴌokas et al. 2008; ᴌokas 2009).

Other interesting aspects that the orbital properties allow us to explore are whether some of the MW satellites were accreted to the MW as part of a group of dwarf galaxies, as is suggested for example by the discovery of possible satellites of the LMC (Koposov et al. 2015a; Bechtol et al. 2015; Martin et al. 2015), and whether the satellites are preferentially distributed on a planar structure (Pawlowski et al. 2012).

These are only a few of the questions that can potentially be addressed by knowledge of the orbital properties of the MW dwarf galaxies. The very first step towards this goal is the determination of the 3D bulk motion of these systems. While determinations of the systemic heliocentric line-of-sight (l.o.s.) velocity are available in the literature for all of the MW classical dwarf spheroidal galaxies and most of the ultra-faint systems discovered so far, astrometric proper motions have became available for the full set of MW classical dwarf spheroidals only very recently and no astrometric proper motion had been available for any of the ultra-faint dwarf galaxies, apart from Segue 1 (Fritz et al. 2018b).

The second release of data from the Gaia mission (Gaia Collaboration 2018a) has started a revolution in this respect. Gaia Collaboration (2018b) demonstrated the power of the Gaia DR2 data for the study of the kinematics of stellar systems around the Milky Way using 75 globular clusters, the Magellanic Clouds, the nine classical MW dwarf spheroidal galaxies, and the Boötes I UFD. In several cases the precisions are exquisite; in others they are comparable to what can be achieved with HST, but with the clear advantage of being in an absolute reference frame.

Gaia DR2 has also opened the door to the determinations of the systemic proper motion of the dozens of other dwarf galaxies surrounding the Milky Way. In this work, we obtain systemic proper motions for most of the systems within 420 kpc of the MW. The paper is structured as follows: in Sect. 2 we describe the data sets used for the analysis and in Sect. 3 the methodology adopted to determine systemic proper motions; in Sect. 4 we present the determination of space velocities and orbital properties of the objects in the sample; we discuss the main results in Sect. 5; and we present the conclusions and summary in Sect. 6.

2. Data

Our sample of galaxies consists of the vast majority of galaxies and galaxy candidates within 420 kpc of the Milky Way. We omitted the Magellanic Clouds because they have well-established motions and their measurable internal motions complicate measurements; Leo T because there are only a handful of stars with spectroscopic data that have matches in Gaia, and these are very faint, and therefore with large proper motion errors; and Bootes III due to its unclear nature as a galaxy or stream (Grillmair 2009). We included objects whose nature as a globular cluster or a dwarf galaxy is still under debate, such as Crater I (Kirby et al. 2015; Voggel et al. 2016).

Our work is mainly based on the Gaia (Gaia Collaboration 2016) data release 2 (GDR2; Gaia Collaboration 2018a). Furthermore, we use literature spectroscopic data sets of dwarf galaxies to ease identification of members (see Table 1).

Because we rely on spectroscopic data, some satellites are omitted from our sample; for instance, we drop Pegasus III because of the lack of Gaia DR2 matches for the spectroscopic members in Kim et al. (2016). Some other systems are so faint and/or distant (e.g., Indus II) that calculations show that no member stars are expected above the Gaia DR2 magnitude limit. Others, like Sagittarius II and Pictor II, are bright enough in principle, but no spectroscopic data set has been published yet.

3. Proper motions

3.1. Method

As a first step, we performed a nonrestrictive selection based on the information from spectroscopic measurements, mainly aimed at excluding obvious contaminants. This was done by retaining stars with at least 40% probability of being members to a given galaxy according to the spectroscopic analysis. Except for a few cases, membership probabilities (binary or continuous probability) are available from the literature. For those cases where probabilities were not provided in the source paper, we calculated them ourselves (see notes in Table 1) using the heliocentric velocities, ν_LOS, and when available also gravity indicators; essentially we assigned probabilities based on the likelihood of the star being a giant and having a velocity similar to the systemic velocity of the galaxy (see Appendix A for details on the target selection).

This first, broad selection of probable members is later refined by including information from Gaia DR2. We first checked whether stars with spectroscopic information have a match within 1″ in Gaia DR2 and also have Gaia DR2 kinematic information. We then required that the stars should have a parallax deviating less than 2σ¹ from the value expected given the distance to the galaxy, and tangential velocities (in right ascension and declination direction) less than twice the corresponding measurement error from an approximation of the escape speed at the distance of the object. Our estimate of the local escape speed corresponds approximately to that expected for a Navarro, Frenk and White (NFW) dark matter (DM) halo (Navarro et al. 1996) with virial mass 1.2 × 10¹² M_⊙. This is a generous selection in which objects with fairly extreme velocities such as Leo I, Boötes III, and the Magellanic Clouds are required to be bound to the MW. After converting the escape speed to mas yr⁻¹ at the distance of the object, we also added the expected “proper motion” due to the reflex motion of the Sun given the position and distance of the satellite. The escape speed criteria was applied separately in the two dimensions (right ascension and declination). All these criteria are relatively inclusive, but because we used several of them the sample is still expected to be rather clean.

Finally, we applied an outlier rejection criteria, excluding those stars whose proper motion deviates by more than 3σ in at least one proper motion dimension, where in this case σ is the measurement error in the individual proper motion.

Overall not many stars were excluded from the initial samples (see Table 1), especially for the fainter galaxies. We also explored changing our spectroscopic probability cut. Even when we used a spectroscopic membership probability > 75%, the change in the resulting values of the systemic proper motion was at most 0.5σ, but typically much less. We list the few exceptions in the appendix. We note that for most of the faint galaxies the membership probability provided in the source papers is binary, and thus the same stars are always selected. A few objects have only two likely members with kinematic information in Gaia DR2 and passing our criteria (see Table 1); therefore, the final step of outlier rejection cannot be applied on these samples containing only two stars. We inspected the motions of these cases with greater care; most of them agree well with each other, and we use them in the following without distinction (see appendix for details). Especially for objects for which only a handful of spectroscopic members are available, we also double checked that Gaia DR2 stars within one half-light radius from the center of these systems, but without spectroscopic measurements, clump close to the member stars in proper motion space. For simplicity we did not add these others stars to our sample.

After we had selected which stars are members, we calculated the average proper motion of the object in right ascension (μ_α^*) and declination (μ_δ) and the average correlation coefficient between μ_α* and μ_δ by taking an error weighted average of all the member stars. These values are given in Table 2. In this table we also define the abbreviations for satellite names used in this paper. For the conversion to velocities and for our analysis we also add a systematic error to the proper motions errors in both dimensions, which we assume to be uncorrelated between α and δ. For this error, we use a value of 0.035 mas yr⁻¹ for galaxies which cover at least 0.2° (see Gaia Collaboration 2018b). For smaller galaxies we use an error of 0.066 mas yr⁻¹ (Table 4 of Lindegren et al. 2018) for galaxies of zero size and we interpolate linearly between these two cases. We list the systematic errors adopted in Table 2. Even though the smallest galaxies tend to have the largest systematic errors, they are always smaller than the statistical errors. On the other hand, for most of the classical dwarf galaxies, the systematic error dominates.

The obtained proper motions are of high quality as shown by the comparison with independent measurements (see Sect. 3.2) and by the reduced χ² distribution between individual proper motion measurements and sample averages, which follows what is expected for samples of the sizes considered here. We note that while in the inner halo there are still many halo stars in streams and in the smooth halo component, we expect the contamination from halo stars to be much lower at large Galactocentric distances. Thus, the motions of distant satellites will be more reliable in that respect.

3.2. Comparison with other measurements

Gaia DR2 proper motions are known to have systematic uncertainties on the order of 0.035 − 0.066 mas yr⁻¹ (Gaia Collaboration 2018b; Lindegren et al. 2018). Since some of the closer galaxies in our sample have measurement uncertainties smaller or comparable to this, it is important to compare our Gaia DR2 proper motions with those measured using independent instruments and methods. Gaia Collaboration (2018b) compared their Gaia DR2 proper motion measurements of nine classical dSphs with those reported in the literature. They found that in general, the Gaia DR2 measurements are consistent with previous ones (based on ground-based and HST data) at the 2σ level. Moreover, when compared to only HST-based measurements the agreement becomes even better, especially when systematic uncertainties are considered. Our proper motions presented in Table 2 are consistent with the results by Gaia Collaboration (2018b) within the quoted 1σ uncertainties for six out of the nine classical dSphs. The exceptions are large galaxies on the sky, where systematic errors and real physical differences might be more important (see Appendix A for details). Thus, our results for classical dSphs are generally in good agreement with Gaia Collaboration (2018b) and we are able to make the same assessments as in their paper.

For the UFDs, the only object for which a comparison to non-Gaia-based proper motions can be made is Segue 1; however, due to a clearly lower precision of this measurement, it is not so useful to judge the goodness of our results (Appendix A). For the other UFDs, there are other Gaia-based measurements available (Simon 2018; Massari & Helmi 2018; Kallivayalil et al. 2018; Pace & Li 2018), which can be used to check the results against possible sources of errors such as selection of members stars, but not for systematic errors. We make a detailed comparison with these works in Appendix A. In most cases we agree better than 1σ, but such a good agreement is to be expected because some of these determinations are not truly independent since they use the same spectroscopic data sets for member selection and as a starting guess. A difference larger than 2.5σ occurs for Crater II and Tucana III compared to Kallivayalil et al. (2018), which is possibly due to a non-optimal treatment of the background contamination. Our measurement disagrees with both of the other two existing measurements for Segue 2 (which also disagree with each other). This galaxy is a difficult case due to its faintness and large foreground contamination, which prevents a safe identification of true members.

4. Velocities and orbital parameters

In order to convert the measured proper motion into tangential velocity in a heliocentric reference frame, we adopted the distance modulus listed in Table 1, giving priority to distance estimates from studies of variable stars, for those cases where multiple works based on different estimators were available (e.g., mean magnitude of the horizontal branch or tip of the red giant branch). As error in the distance modulus, including that from the original sources (see Table 1), we added a floor of 0.1 mag in quadrature as a safeguard against underestimated systematic errors. The l.o.s. and tangential heliocentric velocities are then transformed into velocities in a Cartesian heliocentric (and then Galactocentric) reference system. From these values we then calculate orbital poles.

Table 2 lists the 3D Galactocentric velocities, ν_3D, and the Galactocentric radial and tangential velocity components, ν_rad and ν_tan; the values were derived assuming a 8.2 ± 0.1 kpc distance of the Sun from the Milky Way center and a velocity of the Sun with respect to the Milky Way center of (11 ± 1, 248 ± 3, 7.3 ± 0.5) km s⁻¹ in a reference frame in which the x-axis is positive pointing inwards from the Sun to the Galactic center, the y-axis pointing in the direction of Galactic rotation, and the z-axis pointing towards the north Galactic cap (Bland-Hawthorn & Gerhard 2016).

We derive the errors of derived properties from Monte Carlo (MC) simulations. Typically, they are normal “forward” MC simulations, i.e., where the simulated quantities are obtained by random extraction from Gaussian distributions centered on the observed value and with dispersion given by the error in the observed quantities. However, for positive defined quantities such as V_tan and V_3D this produces a bias when the errors are close to being as large as or larger than the expected values, with the median of the distribution of values from the MC simulations not coinciding with the observed “error-free” value. We take care of this aspect by running “backward” MC simulations (see Appendix B for details). We show the total velocities of all galaxies in Fig. 1.

Fig. 1. Total velocity of all galaxies in the sample (points with error bars). The curves show the escape velocity for the two potentials used in this paper, i.e., MWPotential14 with a NFW halo of virial mass 0.8 × 1012 M⊙ (black line) and a more massive variant with virial mass 1.6 × 1012 M⊙ (red line).

We then use the publicly available code galpy to integrate orbits in a Galactic potential. We adopt the standard MWPotential14 model, which has three components: a spherical bulge, a disk, and a NFW halo (see Bovy et al. 2015 for details). For the MW dark matter halo mass we explore two cases: a light halo with virial mass 0.8 × 10¹² M_⊙ (Bovy 2015) and a heavier model, identical to MWPotential14 but for the MW DM halo virial mass, which is 1.6 × 10¹² M_⊙ (Fritz et al. 2017). In Table 3 we list the apocenter r_a, pericenter r_p, and eccentricity values derived within these two potentials. The eccentricity is calculated as ecc = (r_apo − r_peri)/(r_po + r_peri) within galpy. The distance of the Sun from the Galactic center and the Sun motion are slightly different than the values used above, but in previous works it was found that this only caused minor changes to the orbital properties; for our sample, the observational sources of errors are dominant in nearly all cases. We use MC realizations of the orbit integrations to estimate the errors on the orbital parameters (see Appendix B.2 for details). Also, because our calculations do not include the effect of the gravitational potential of M 31, the results for systems with very large apocenters should be taken with a grain of salt and should be seen as an indication that the given object can be on its first infall to the Milky Way.

Figure 2 shows the orbital properties of the objects in the sample, quantified as the pericenter versus apocenter (for the heavy and light NFW halo in the top and bottom panel, respectively). We only plot objects with an error in 3D velocity of less than 100 km s⁻¹; this is less than half the escape speed expected for the low-mass MW-halo at the largest distances probed (see Fig. 1) and avoids biases towards large values of apocenter for galaxies with larger errors.

Fig. 2. Pericenter vs. apocenter properties of the galaxies in the sample with 3D velocity errors < 100 km s−1, as indicated by the labels. The top panel plot shows the results for a MW DM halo of virial mass equal to 1.6 × 1012 M⊙, and the lower panel for 0.8 × 1012 M⊙. We also show typical virial radii in both plots. Gru I is omitted in the lower plot because of the very large apocenter.

In Fig. 3 we plot an all-sky view of the orbital poles of the objects in the sample, concentrating on those within 200 kpc, and comparing their location on this plane with the vast polar structure (VPOS) of satellites (Pawlowski et al. 2012).

Fig. 3. All-sky view of orbital poles for the objects in the sample; the circles indicate the median of the 2000 Monte Carlo simulations, while the small points around each object plot the orbital poles from the individual simulations. The magenta circles contain 10% of the sky around the assumed VPOS pole, which is given as an “X” for the co-orbiting direction (orbital sense as most classical satellites, including the LMC and SMC), and a “+” for the opposite normal direction (counter-orbiting). Top panels: objects with Galactocentric distances between 0 and 50 kpc (left) and 50 and 100 kpc (right); bottom panels: objects with Galactocentric distances between 100 and 200 kpc (left) and in the putative Crater-Leo group (right).

5. Results and discussion

5.1. Apocenter distances and orbital energies

In Fig. 1 we show how the newly measured 3D velocities of the satellites as a function of radius compare to the local escape velocity in the two MW potentials. All the galaxies analyzed are compatible with being bound to the MW, in both the light and heavy MW potentials. This includes very distant systems like Phoenix and Eridanus II, although the large error bars could place them away from the locus of bound galaxies. The proper motion errors for both are still too large to change the result based on LOS velocities, that Eri II is likely bound (Li et al. 2017) and that Phoenix I is likely on first infall (Kacharov et al. 2017).

Figure 2 allows us to see the number of galaxies that have orbits that will take them beyond the virial radius of the MW DM halo, i.e., those that have apocenters above the solid lines. These satellites can be considered alongside the backsplash population of satellites that exist beyond the virial radius at z = 0 in simulations, but that were within the virial radius at earlier times (Balogh et al. 2000; Gill et al. 2005; Teyssier et al. 2012; Simpson et al. 2018). Considering only the objects with errors in ν_3D <  100 km s⁻¹ (30 systems), Leo I is presently located beyond the virial radius in both MW potentials, while CVen I and Leo II only in the light MW halo. When considering the values within the 1σ range, only Leo I is found to be certainly “backsplashing” in the heavy halo model, with other possible candidates being CVen I, Leo II, UMa I, Gru I, and Boo II; the last two might fall in the category of systems infalling for the first time, given their negative Galactocentric radial velocities.

In the low-mass halo, 13 objects have apocenters that would still exceed the virial radius even after subtracting the 1σ errors, and are therefore backsplashing satellites in that their orbits will take them beyond the virial radius (5 of these objects have negative ν_rad); several others most likely have apocenters smaller than the virial radius, but could potentially be scattered beyond it due to the measurement uncertainty. This would imply that a large fraction of the observed dwarfs could be backsplashing, while for a high-mass halo the fraction is very much reduced.

The definition of whether or not a satellite is bound in an evolving cosmological context is complex, and the comparison to the curves in Fig. 1 should therefore be interpreted primarily in a relative sense. Boylan-Kolchin et al. (2013) addressed this topic in detail for the case of Leo I. That analysis was based on the proper motion of Leo I measured by HST, which is in good agreement with our Gaia DR2 value. They showed that the low-mass MW is actually ruled out at 95% confidence, since in such halos it is vanishingly rare in cosmological simulations to find subhalos at 273 kpc moving as fast as Leo I. By contrast, they showed that the high-mass Milky Way is the statistically preferred value. Given this result, our Figs. 1 and 2 show that several other satellites (e.g., Car II, Car III, Gru I) have a combination of distance and velocity that will put them at odds with a low-mass MW.

5.2. Pericenter distances and tidal influences

From Fig. 2 it can be seen that the measured proper motion of Tuc III leads to a very internal orbit in both MW potentials, with the object reaching well within 10 kpc and being confined within 50 kpc from the MW center. Tuc III is therefore highly likely to have been subject to strong tidal disturbances, as confirmed by the presence of a stream around it and predicted by Erkal et al. (2018). Nonetheless, the proper motion we measure is different from the values predicted in Erkal et al. (2018), independently of their adopted LMC mass; this could be indicative of a difference between the MW potential assumed in that work and the actual MW potential. We note that Tuc III is one of the few satellites in our sample whose nature as a galaxy is not confirmed according to the spread in the distribution of l.o.s. velocities and metallicities (Simon et al. 2017; Li et al. 2018a); perhaps its possible nature as a stellar cluster could explain why its pericenter is much smaller than those of all the other systems.

The following galaxies have a likely pericenter smaller than 20 kpc in at least one of the MW potentials: Wil 1, Seg 1, Sgr I, Cra II, Dra II, Tri I, and Her I. This is clearly a dangerous region, as empirically highlighted by the obvious state of tidal disruption afflicting Sgr I (Majewski et al. 2003). The other galaxies broadly consist of two groups. The first contains relatively bright galaxies (Cra II, Her I), which also show signs of tidal disruptions (see below). The second group contains very faint satellites which, due to selection biases, can only be detected if currently relatively close to us, and thus have a higher probability of also having a smaller pericenter (Wil 1, Seg 1, Dra II, Tri II). Apart from Dra II (Longeard et al. 2018), for which efficient cleaning of foreground contamination was done through the use of narrow-band filters, no signs of tidal disturbance have been detected yet in the other three systems. This could be due to the intrinsic difficulty in unveiling such signs in extremely faint systems in the presence of contamination, or could also be a consequence of their relative compactness (with half-light radii of at most 21 pc; Carlin et al. 2017). We note that, within the more massive MW potential, the galaxies with observed tidal features have closer pericenter passages with respect to the case of the less massive MW, another argument in favor of a more massive halo. The inhomogeneous nature of galaxies found to have small pericenter distances might also be partially caused by large uncertainties in the derived quantities.

As mentioned in the Introduction, elongations and/or isophote twists suggestive of tidal stretching of the stellar component have been observed in Carina I, Hercules, and Boötes I. For Her I a pericenter as small as 5–10 kpc is within the 1σ bounds for both potentials, thus it supports the nature of the observed elongations being of tidal origin. Our proper motion is within 0.7/1σ of the proper motion prediction of Küpper et al. (2017), who explain the structure of Hercules on an orbit with a pericenter distance of about 5 kpc.

On the other hand Carina I is not expected to have come closer than ∼50 kpc (within 1σ) for the heavy halo, while in the case of the light halo its orbit would be rather external. The observational findings of signs of likely tidal disturbance in Carina can be considered as robust as they have been detected in multiple studies in the literature adopting different methodologies (see, e.g., Muñoz et al. 2006; Battaglia et al. 2012; McMonigal et al. 2014, and references in); their presence might then suggest a preference for a heavy MW DM halo or that dwarf galaxies might experience strong tidal disturbance even on rather external orbits. We deem this latter hypothesis unlikely (see results from N-body simulations, e.g., Muñoz et al. 2008; Peñarrubia et al. 2008, for those simulated objects on similar orbits), but this is certainly an important aspect to be verified since several of the MW dwarf galaxies are on more internal orbits than Carina (Draco and Ursa Minor to mention two) and could see their internal kinematic properties potentially affected. We note though that the hypothesis put forward by Fabrizio et al. (2016) that Carina is the result of an initially disky and rotating dwarf galaxy heavily tidally stirred by the Milky Way would require it to be on an exactly prograde and tighter (125 kpc/25 kpc apocenter/pericenter) orbit than measured here even with the heavy MW potential, and it does not therefore seem likely. Alternatives could contemplate the possibility that the properties of Carina might have been affected by the Magellanic system: e.g., the recent infall of the Magellanic system might have modified Carina’s current orbital properties, making them therefore not representative of its past orbital history; however, Carina and the LMC have very similar orientations of the time-average orbit (Gaia Collaboration 2018b), while torques from the LMC would be expected to be maximal on systems on a perpendicular orbital plane; we speculate that Carina could have experienced a close encounter with the Magellanic Clouds and that perhaps these latter systems are responsible for inducing tidal disturbances. While the perturbative influence of the LMC (Gómez et al. 2015) has not been included in the presented orbit calculations, it will be included in follow-up studies to assess the impact on the orbits of the classical and UFD satellites of the MW (Patel et al., in prep.).

Depending on the potential adopted, the Boötes I case is intermediate between Car I and Her I, with a pericenter ≳20 kpc but for which smaller values cannot be excluded (see also Gaia Collaboration 2018b; Simon 2018). Thus, the observed elongated and structured spatial distribution of its stars could still be compatible with tidal features caused by the MW.

Crater II is a remarkable galaxy; it is larger than Fornax, but only as luminous as CVen I (Torrealba et al. 2016a) and has a very cold internal kinematics when compared to galaxies of similar size (Caldwell et al. 2017). Its peculiar properties have raised interest in the community. Its low l.o.s. velocity dispersion appears compatible with the predictions given by MOND based only on the characteristics of the Cra II stellar component (McGaugh 2016), at least under the assumption of dynamical equilibrium. On the other hand, its internal kinematical and structural properties could also be consistent with Crater II having been embedded in a typical dark matter halo expected in ΛCDM for dwarf galaxies, which has undergone heavy tidal stripping (see, e.g., Fattahi et al. 2018; Sanders et al. 2018).

In Fig. 4 (top) we plot Cra II distance from the Galactic center as a function of time: in both MW potentials, the best Cra II orbit is radial, with an eccentricity of about 0.7, and can reach as close as 2 kpc from the MW center. That implies that the Hill radius r_Hill = (m_CraterII/(3 * m_MW))^1/3 * r_{peri, CraterII} (where m_MW is the MW mass within the pericenter distance of Crater II, r_{peri, CraterII}) is far smaller than Cra II half-light radius (see Fig. 4, bottom panel), therefore confirming the prediction of Fattahi et al. (2018) and Sanders et al. (2018) that it is in process of tidal disruption. Our proper motion is in good agreement (0.3/1.3σ) with the models by Sanders et al. (2018).

Fig. 4. Properties of Crater II. Top: orbital integration for Cra II in the MWpotential14 for the standard light halo (blue line) and the heavy DM halo (red line) plotted as the distance from the MW center as a function of time. Bottom: hill radius of a galaxy of the mass of Cra II (Caldwell et al. 2017) in two different MW potentials compared to the size and apocenter of Cra II.

As stated above, several other systems have internal orbits similar to the orbit of Crater II and therefore the potential impact of tidal disturbances needs to be understood to properly interpret these systems’ structural and internal kinematic properties.

5.3. Missing satellite problem

Figure 5 shows histograms of the ratio f = (d_GC − r_peri)/(r_apo − r_peri): a galaxy with f close to zero is near pericenter, while a value close to one indicates that the galaxy is near apocenter. Both histograms have a weak peak at small value of f. This is more pronounced for the low-mass MW model; based on the arguments in the preceding sections, we consider this less plausible. We also show the variant where only galaxies with ν_3D <  100 are plotted since we do not expect their orbital properties to suffer from significant biases. When galaxies with large errors are excluded, the histogram for the higher mass halo is close to flat. However, basic dynamics dictates that, within their orbits, galaxies spend most of their time near apocenter where the velocity is lower. If the orbits were circular the argument would not be valid, but the median eccentricity is 0.53 for the high-mass MW halo, when the galaxies with larger errors are excluded.

Fig. 5. Histograms of the ratio f = (dGC − rperi)/(rapo − rperi) for high-mass and low-mass MW models (upper and lower panels, respectively). The galaxy on the left of each histogram is near pericenter, while the galaxy on the right of each histogram is near apocenter.

Also, the number of galaxies at small pericenters (≲20 kpc) is reduced through tidal destruction. If we had a complete sample of MW dwarfs, then the histograms would have to be increasing towards high f. By contrast, even for the high-mass MW, the observed histogram is flat at best. We could argue that the peak towards lower values of f could be caused by a group of former satellites of a larger galaxy (e.g., the LMC) having infallen; however, the number of objects that we find as possibly having been associated with each other is not sufficient to explain the feature (see Sect. 5.5). The corollary is that there must be a population of ultra-faint dwarf galaxies that are currently at apocenter, especially beyond ∼100 kpc, that have yet to be discovered. This is relevant to the so-called missing satellite problem (see review in Kravtsov 2009) as it affects the comparison of observed dwarf galaxy counts to subhalos found in hierarchical galaxy formation scenarios. The implications of this statement will be discussed in Patel et al. (in prep.).

5.4. Orbital poles and planar alignments

Figure 3 shows the distribution of orbital poles for three distance bins (0–50, 50–100, and 100–200 kpc), and for the proposed members of the Crater-Leo Group (not included in the 100–200 kpc plot). The uncertainty in the direction of orbital poles is illustrated with point clouds based on the 2000 Monte Carlo simulations of the measurement uncertainties. The magenta cross and plus sign give the assumed normal direction to the VPOS, as used to predict the proper motions of satellites in Pawlowski & Kroupa (2013). The magenta circles of opening angle θ_inVPOS = 36.87° denote areas of 10% of the sphere around these directions. We consider orbital poles that lie within this region to orbit along the VPOS.

Based on the assumption that the satellite galaxies orbit along the VPOS, Pawlowski & Kroupa (2013) and Pawlowski et al. (2015) predicted orbital poles for these systems. We can now test how well these predictions are met. The predicted orbital pole direction is the direction along the great circle perpendicular to the satellite (as seen from the center of the Milky Way) that minimizes the angle to the assumed VPOS normal, which points to Galactic coordinates (l, b)=(169.3 ° , − 2.8 ° ). The corresponding minimum angle between the VPOS normal and the pole of a satellite galaxy is θ^pred. The angle between the orbital pole based on our measured proper motions and the VPOS normal is θ^obs. The ratio between these gives a measure of how well an observed orbital pole agrees with its predicted direction. Figure 6 plots this ratio against a measure of the uncertainty in observed orbital pole direction, Δ_pole. The latter is defined as the angle from the most-likely measured orbital pole of a satellite which contains 68% of the orbital poles sampled from its measurement uncertainties. The ratio of observed to measured angle from the VPOS is typically high for objects with poorly constrained orbital pole directions (Δ_pole >  30°), but for those with smaller uncertainties the majority of objects (17 of 24) agree to within a factor of two (dashed magenta line) with the predicted angle.

Fig. 6. Alignment with the VPOS, measured as the ratio between the measured (θobs) and the predicted (θpred) angular offset of an orbital pole from the VPOS, plotted against the uncertainty (Δpole) in the direction of the orbital pole. Blue upward triangles are co-orbiting relative to the majority of classical satellites including the LMC; red downward triangles are counter-orbiting. Symbol size increases with Galactocentric distance.

Table 4 provides an overview of the predicted (θ^pred) and the measured (θ^obs) alignments of satellite orbital poles. The observed angle is given a negative sign if the measured orbital pole is counter-orbiting relative to the orbital direction of the majority of the classical satellite galaxies associated with the VPOS. For each satellite, we calculate the fraction p_(inVPOS) of Monte Carlo sample orbital pole directions which falls to within θ_inVPOS of the VPOS normal. For an orbital pole that is misaligned with the VPOS, this indicates whether there is nevertheless a chance that the pole might be aligned. Also given in the table is a measure of how strong the constraints of each satellite’s pole are on its alignment with the VPOS direction. For the latter, we assume that the satellite’s intrinsic orbital pole aligns perfectly with its predicted direction. We then vary its orbital pole direction by sampling 2000 times from the measured uncertainties in the orbital pole direction of this satellite. Since a satellite’s orbital pole cannot be better aligned than the predicted direction, any offset from this direction results in an increase in the angle with the VPOS normal. For each realization, the angle to the VPOS normal is calculated. If it is larger than θ_inVPOS, we count this realization as not aligned with the VPOS, even though we know that its intrinsic pole is perfectly aligned. The fraction of realizations counted this way gives an estimate of the probability p_{> VPOS} of falsely finding this satellite’s orbital pole to be misaligned with the VPOS. It is compiled in Col. 5 of the table. We also count how often the angle to the VPOS exceeds θ^obs[°], the angle between the most-likely measured orbital pole and the VPOS normal. This gives an estimate of the probability p_{> obs} of measuring an intrinsically well-aligned orbital pole as far away from the VPOS normal as observed. Since this method assumes intrinsically perfect alignments, the resulting probabilities should be seen as lower limits.

Overall, for 12 of the satellites the chance p_inVPOS of aligning with the VPOS to better than θ_inVPOS given our proper motion uncertainties is lower than 5%. For the remaining 27 satellites, an alignment of the orbital pole with the VPOS is either found, or cannot be rejected with high confidence. Six of the 39 objects cannot align with the VPOS because their predicted angle θ^pred already exceeds θ_inVPOS: their spatial positions alone already place them outside the VPOS plane orientation. A well-known example is Sagittarius, which has an orbit almost perpendicular to the VPOS. The other satellites in this category are Hercules I, Segue 2, Triangulum II, Ursa Major II, and Wilman 1. Five of these satellites counter-orbit relative to the VPOS, for a counter-orbiting fraction of $f_{\mathrm{counter}}=\frac{5}{6}=0.83$. Of the remaining 33 satellites for which an alignment is feasible, 17 have median orbital poles that align to within θ_inVPOS with the VPOS normal. The majority of these are well constrained to align with the VPOS because most of their Monte Carlo sampled orbital poles also show an alignment (p_inVPOS >  0.5). The three least certain alignments are CVen I, Leo V, and Dra II, for which only 38%, 40%, and 48% of realizations fall to within θ_inVPOS of the VPOS normal, respectively. Of these 17 satellites with aligned orbital poles, only 1 has p_{> VPOS} ≥ 0.5 (Leo V), 5 have 0.4 ≥ p_{> VPOS} ≥ 0.17, and the others have p_{> VPOS} <  0.1. This indicates that given the respective measurement uncertainties, the measured orbital poles are expected to be found within θ_inVPOS of the VPOS if the satellites orbit along the VPOS. It is also interesting to note that of these 17 satellites, 6 counter-orbit: Aqu II, CVen II, Leo V, Seg 1, Tuc III, and Scu I (for which this was known previously). The remaining 11 satellites co-orbit, as do the LMC and SMC, so the counter-orbiting fraction is $f_{\mathrm{counter}}=\frac{6}{17+2}=0.32$.

The remaining 16 satellites have median orbital poles that do not align with the VPOS, 9 of which have orbital pole directions that are too weakly constrained to be conclusive. These satellites are Cra I, Eri II, Gru I, Hya II, Leo I, Leo II, Leo IV, Phx I, and Pis II. When the orbital pole directions are Monte Carlo sampled from the measured uncertainties, there is about a one in three chance that the pole is aligned with the VPOS (0.13 <  p_inVPOS <  0.43), thus consistent with aligning with the VPOS within their uncertainties. Even if their poles were intrinsically as well aligned as geometrically possible, due to the considerable uncertainty in their orbital pole directions their poles would be expected to be found outside of θ_inVPOS with probabilities of p_{> VPOS} = 37%–60%. With p_(>obs) = 10%–51%, it is expected that the poles are as far as or farther from the VPOS as observed, except for Leo I and Gru I for which these probabilities are only 1.1% and 0.4%, respectively. Three of the seven are more likely counter-orbiting, the remaining five more likely co-orbiting relative to the VPOS, so the counter-orbiting fraction is $f_{\mathrm{counter}}=\frac{3}{7}=0.43$.

Of the remaining seven satellites that have orbital poles that are misaligned with the VPOS normal, one is consistent with aligning to within θ_inVPOS at a 6 % level (Boo I). However, it is basically impossible that this object aligns as well as geometrically possible with the VPOS because p_{> VPOS} ≈ 0. This means that if this satellite orbits in the VPOS, then its pole cannot be quite as close to the VPOS normal as is geometrically possible. The satellite is most likely counter-orbiting.

Thus, only six of the satellites have orbital poles that are firmly and conclusively misaligned with the VPOS even though they could have had aligned poles: e Boo II, Car II, CBer I, Sxt I, Tuc II, and UMa I. They all have well-constrained orbital pole directions p_inVPOS ≤ 0.02 and p_{> obs} = 0, thus are not orbiting along the plane defined by the VPOS. Three of the six are most likely counter-orbiting (f_counter = 0.5). For a more in-depth analysis in regard to the VPOS, we refer the reader to Pawlowski et al. (in prep.).

5.5. Group infall

The orbital pole of the LMC is similar to the direction of the VPOS normal. Therefore, our proper motion measurements suggest that several of the dwarf galaxies in the sample that have orbital poles falling within 10% of the sky around the assumed VPOS normal in the co-rotating assumption might be associated with the LMC. The association of MW dwarf galaxies with the Magellanic system is the subject of other works (Kallivayalil et al. 2018; Fritz et al. 2018a).

Mainly on the basis of the objects’ position on the sky and heliocentric distance, Torrealba et al. (2016a) argued that Crater I, Crater II, Leo II, Leo IV, and Leo V might have once formed part of a group accreted by the MW. Even though the errors on the orbital pole determinations are large, our analysis suggests that such a prior physical association of all five objects together is unlikely. Given the large measurement uncertainties, the poles of Leo IV and V are rather unconstrained, but compatible with each other. An association of Cra I and Cra II appears likely, since several of the MC realizations bring the Cra I orbital pole in agreement to the region occupied by Cra II pole, and the two systems share very similar orbital properties (or compatible within 1σ).

In contrast the poles of Cra II and Leo II, which are better constrained, do not overlap well. The preferred pericenter of these two objects differs by about 30 kpc, but the values agree within 1σ; therefore, it cannot be used as a discriminator. Also, the orbital properties of Cra I are compatible with those of Leo II.

Thus, we cannot exclude the possibility that there are two groups, one including Cra I and Cra II, and the other with Leo IV and Leo V, with Leo II being possibly part of one group or the other; we deem it unlikely that all five objects came together as one group. However, since the HST proper motion measurements places the orbital pole of Leo II in the VPOS direction (Sohn et al. 2013; Piatek et al. 2016; Pawlowski et al. 2017), it might be possible that none of the measurements of Leo II is accurate enough and that further studies with either HST or Gaia are necessary.

6. Summary and conclusions

We derive systemic proper motions for all dwarf galaxies or galaxy candidates within 420 kpc using Gaia DR2 proper motions for which literature spectroscopic members are available. Our proper motion determinations are in very good agreement (usually within 1σ) with the determinations by Gaia Collaboration (2018b), Simon (2018), Massari & Helmi (2018), and Pace & Li (2018).

We derive the implied Galactocentric velocities, and calculate orbits in canonical MW halo potentials of low (0.8 × 10¹² M_⊙) and high mass (1.6 × 10¹² M_⊙). Comparison of the distributions of orbital apocenters and 3D velocities to the halo virial radius and escape velocity, respectively, suggests that the satellite kinematics are best explained in the high-mass halo. Relevant to the missing satellite problem, the fact that fewer galaxies are observed to be near apocenter than near pericenter implies that there must be a population of distant dwarf galaxies yet to be discovered.

Several satellites have likely pericenter distances reaching within 20 kpc of the Milky Way centers, and are therefore candidates for having suffered strong tidal disturbance (Tuc III, Sgr I, Her I, Dra II, Cra II, Tri II, Seg 1, and Wil 1). Among these satellites, the orbital properties of Tuc III, Sgr I, Her I, Dra II, and the “feeble giant” Crater II are in line with the detections of tidal features in these objects and the predictions of models. In contrast, no tidal features have been detected yet in Tri II, Seg 1, and Wil 1. This difference could be partly caused by the still relevant errors in the derived pericenter values and/or by the intrinsic difficulty in detecting tidal features, especially in very faint systems. Our analysis suggests that some classical MW dSphs, such as Draco I and UMi, also have internal orbits with pericenter distances that bring them dangerously close to the internal regions of our Galaxy, at risk of being tidally affected. It appears then crucial to address in more detail the impact that tidal disturbances might have had in the structural and internal kinematic properties of these galaxies.

Of the 23 satellites for which we can draw conclusions, 17 are orbiting along the plane of satellites (VPOS, not counting the well-aligned LMC and SMC) and 6 are not. These findings suggest that a majority of the MW satellites for which we have measured proper motions orbits along the VPOS, but that not all satellites participate in coherent motion along this structure. This is in line with an analysis based solely on the spatial distribution of satellite galaxies, which found that up to half of the MW satellites might be drawn from an isotropic distribution in addition to satellites drawn from a planar distribution (Pawlowski 2016), and with the satellite plane around Andromeda which consists of about half of the M 31 satellites (Ibata et al. 2013).

The distribution of orbital poles does not appear to confirm the hypothesis that Crater 1, Crater II, Leo II, Leo IV, and Leo V were all accreted by the MW as part of the same galaxy group; although the errors are large for most and only a common infall of Crater II and Leo II is made clearly less likely by our measurement.

Finally, we note that we just use part of the power of Gaia DR2 for the determination of systemic proper motions of dwarf galaxies and this has already led to constraining proper motions for dozens of galaxies. We expect that in some cases adding stars without existing spectroscopic measurements should improve the precision. Since the precision in proper motion determinations grows with the 1.5 power of the time-baseline (and also systematics, when they are not based on a moving reference frame) we expect that the proper motions of Gaia should be a factor of 4.5 more accurate after the nominal mission and possibly a factor of 12 after the extended mission. This would enable us to measure systemic proper motions for essentially all galaxies in the sample with the best precision now possible for very few galaxies. For example, even the motion of Phoenix I could be measured to about 19 km s⁻¹ precision, similar to the third best ² measurement of HST to date.

Acknowledgments

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. G.B. gratefully acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the Ramon y Cajal Programme (RYC-2012-11537) and the grant AYA2014-56795-P. M.S.P. acknowledges that support for this work was provided by NASA through Hubble Fellowship Grant #HST-HF2-51379.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. N.K. is supported by the NSF CAREER award 1455260. This project is also part of the High-resolution Space Telescope PROper MOtion (HSTPROMO) Collaboration, a set of projects aimed at improving our dynamical understanding of stars, clusters, and galaxies in the nearby Universe through measurement and interpretation of proper motions from HST, Gaia, and other space observatories. We thank the collaboration members for sharing their ideas and software. We thank Josh Simon, Nicolas Martin, and Matthew Walker for providing the catalogs of Simon & Geha (2007), Martin et al. (2016), and an improved catalog of Walker et al. (2016), respectively. We thank Matteo Monelli for guidance on the bibliography for distance estimates from variable stars.

References

Appendix A: Details on galaxies

In this appendix we first give details on aspects concerning all or several of the objects in the sample, and then provide specific details for each object.

• For the systematic error in proper motion, we need to estimate the size of each galaxy or candidate galaxy in the sample. We determine this by taking the average extent of member stars in RA (corrected by cos(Dec)) and Dec and averaging between these two values. For example, for Aqu II this leads to a size of 0.019° and thus a systematic error of 0.063 mas yr⁻¹.

• For spectroscopic catalogs not containing membership information, we adopt different selection criteria, depending on whether only l.o.s. velocities are available from the literature or also other additional information that allows us to classify a star as a giant, such as the star’s surface gravity, logg, or the gravity-sensitive criterion in Battaglia & Starkenburg (2012), based on the equivalent width (EW) of the NIR Mg I line as a function of the EW of the near-IR Ca II triplet. The gravity or gravity-sensitive criteria are only applied to systems for which the spectroscopic members are indeed expected to be giant stars.

When only l.o.s. velocities are used, we assign probabilities (p) with the equation

$$\begin{array}{c}\hfill p=exp(-1/a\ast {({\nu }_{\text{LOS}}-{\nu }_{\text{LOS,sys}})}^2),\end{array}$$

while when logg estimates are also available we use

$$\begin{array}{c}\hfill p=(tanh([c-log(g)]/b)+1)/2\ast exp(-1/a\ast {({\nu }_{\text{LOS}}-{\nu }_{\text{LOS, sys}})}^2),\end{array}$$(A.1)

where ν_LOS is the heliocentric l.o.s. velocity of a given star; ν_LOS, sys the systemic velocity in the same system; and a, b, and c are normalization constants that depend on the system and are chosen to yield a probability of ∼0.5 at heliocentric velocities about 3x the internal l.o.s. velocity dispersion of the system, as given in the source papers.

For Fornax and Sculptor, when considering the Battaglia & Starkenburg (2012) catalogs, we adopt their same selection of likely giant stars.

A.1. Aquarius II

For Aquarius II just the two brightest members (according to the binary classification) of Torrealba et al. (2016b) have matches in GDR2. However, we consider its systemic motion reliable because the star density in the stellar halo is low at the Aquarius distance (105 kpc) and because three stars with similar colors without spectroscopy have a consistent proper motion (see also Kallivayalil et al. 2018).

A.2. Boötes I

We use two spectroscopic catalogs for this galaxy, Martin et al. (2007) and Koposov et al. (2011). The former contains binary membership classification and we also adopt from that study the values of the heliocentric l.o.s. systemic velocity (ν_LOS, sys = 99 km s⁻¹) and velocity dispersion (5 km s⁻¹); no membership classification is given in the latter.

For stars that were observed in both works, we adopt the binary classification of Martin et al. (2007). To those that are only in the catalog of Koposov et al. (2011), we assign probabilities on the basis of the l.o.s. velocity. Our selection is somewhat different from that of Simon (2018) but we agree within 1σ with the systemic proper motion. We also agree within 1σ with the PM determination by Gaia Collaboration (2018b), who used a selection based on photometry.

A.3. Boötes II

We find matches in GDR2 for four bright members (according to the binary classification) in Koch et al. (2009). We consider Boötes II systemic PM reliable since some other stars whose colors are compatible with Boo II share the same motion. Our determination agrees within 1σ with that by Simon (2018).

A.4. Canes Venatici I

We use the catalogs of Martin et al. (2007) and Simon & Geha (2007) and the binary probabilities therein. We give the classifications of Martin et al. (2007) higher priority. Since 57 stars pass all our cuts, the derived motion is robust.

A.5. Canes Venatici II

We use the catalog of Simon & Geha (2007). From its binary probabilities 11 members remain after all cuts.

A.6. Carina I

We use the catalog of Walker et al. (2009b) with continuous probabilities. Our proper motion agrees within 1σ with that of Gaia Collaboration (2018b). For both works the systematic error dominates over the statistical uncertainty.

A.7. Carina II

We use the catalog of Li et al. (2018b), which has binary probabilities. The 18 member stars we select form a clear clump in proper motion space supported by stars without spectroscopy. Our motion agrees within 1σ with Simon (2018), Kallivayalil et al. (2018), and Massari & Helmi (2018).

A.8. Carina III

Again we use the catalog of Li et al. (2018b), which has binary probabilities. The four member stars we select form a clear clump in proper motion space supported by one bright star without spectroscopy. Our motion agrees within 1σ with Simon (2018) and Kallivayalil et al. (2018).

A.9. Coma Berenices I

We use the catalog of Simon & Geha (2007) with its binary membership classifications. After all cuts this leads to 17 members in our final selection. They form a clear clump in proper motion space supported by a few stars without spectroscopic information. Our proper motion agrees within 1σ with the motion of Simon (2018).

A.10. Crater I

For this system we give higher priority to the continuous probabilities by Voggel et al. (2016). If that is not available, we use the binary classification of Kirby et al. (2015). Ten members stars remain after the membership selection. Due to the relative high surface brightness of Crater I, pollution by contaminant stars in not an issue. Because the available spectroscopy is deep and almost complete there are no additional member candidates with Gaia DR2 information only.

A.11. Crater II

In this case, probabilities are not included in the tables of Caldwell et al. (2017); therefore, we calculate them using l.o.s. velocities (with a systemic velocity of 87.5 km s⁻¹) and logg. For log(g) our selection is guided by Fig. 3 of Caldwell et al. (2017). With our selection, 59 stars have a probability greater than 0.4, which is slightly less than in Caldwell et al. (2017) where 62 stars have a probability greater than 0.5, but is on the safe, conservative side. If we were to change the probability cut to 75%, this would yield a 0.56σ difference in proper motion for the declination component. This is 0.019 mas yr⁻¹, less than the systematic error; therefore, we consider our motion as robust.

Recently, Kallivayalil et al. (2018) added 59 stars to our spectroscopic sample, selecting them from photometry, and obtained a 4σ different proper motion. Most of the added stars are rather faint, thus it is not clear that adding them improves the overall accuracy, since it is difficult to be sure that these faint stars are truly members, especially since Crater II has a very low surface brightness. It could also be that assignment of binary membership is not sufficient for photometric members, and that a probabilistic treatment like that in Pace & Li (2018) is advisable. Also, there seems to be a problem for Tuc II with the motion of Kallivayalil et al. (2018; see Sect. A.34). Therefore, we prefer the proper motion derived here.

A.12. Draco I

We use the catalog of Walker et al. (2015). This contains probabilities and we calculate them using the l.o.s. velocities and star’s gravity. Our motion agrees within 1σ with Gaia Collaboration (2018b).

A.13. Draco II

We use the catalog of Martin et al. (2016) and find six stars matching all our criteria. The derived PM agrees very well with that by Simon (2018). There are a few faint stars without spectra whose motion agrees with the spectroscopic selected members (Kallivayalil et al. 2018).

A.14. Eridanus II

We use the catalog by Li et al. (2017) with binary probabilities, 12 of which are bright enough to have Gaia DR2 kinematic measurements. Due to the galaxy faintness and large distance, the errors on the PM are large.

A.15. Fornax (I)

We use the catalogs of Walker et al. (2009b) and Battaglia & Starkenburg (2012). For both data sets we use the LOS velocities to select members for the Battaglia & Starkenburg (2012) data set also the gravity indicators. Due to the brightness of that galaxy, our proper motion does not depend on the membership probability cut adopted. Our proper motion differs in Dec by about 2σ (our statistical error), or 0.012 mas yr⁻¹ from Gaia Collaboration (2018b). Since the spatial coverage of the data we use differs from that in the Gaia collaboration analysis, this discrepancy could be caused by spatial variations in the systematic error or by real physical differences over the angular area covered by this rather large galaxy.

A.16. Grus I

For Grus I we use the catalog by Walker et al. (2016). We note that the probabilities in the long electronic Table 1 of Walker et al. (2016) are incorrect for Grus I. We use a corrected version kindly provided by the lead author. The probabilities listed in Table 1 in the PDF version of that article are correct and contain all members except for one rather faint star. Grus I is one of the few satellites whose systemic PM is somewhat sensitive to which stars are included or excluded. The two brightest stars of its five likely members have p = 0.69 (Gru1-032) and p = 0.68 (Gru1-038), i.e., relatively low membership probabilities. Excluding these two stars changes the velocity by 0.7 and 0.9σ, respectively, in the two dimensions. We decided to include these stars since both their PMs are within our proper motion selection box even without considering their error bars. Our systemic PM agrees with that in Kallivayalil et al. (2018), which is unsurprising since we mostly use the same stars. A motion also agrees with the purely photometric proper motion determination of Pace & Li (2018).

A.17. Hercules (I)

In the case of Hercules we assign the highest priority to the work of Simon & Geha (2007). We complement the analysis with the stars in Adén et al. (2009), which were all classified as members in that study even though not all of them have spectroscopic information, e.g., some candidate horizontal branch. Since the stars without spectroscopy are faint, excluding them causes significant changes in the average velocity, given the error bars. We therefore include them in our final estimate.

A.18. Horologium I

For Horologium I we use the catalog of Koposov et al. (2015b), which contains binary membership classification. Four members remain after all our cuts. These relatively bright stars form a tight group in proper motion space. Our proper motion agrees very well with that in Simon (2018), Kallivayalil et al. (2018), and Pace & Li (2018).

A.19. Hydra II

We use the catalog of Kirby et al. (2015), which contains binary membership classification. After our cuts we have six members, whose motions agree well with each other. Our systemic PM also agrees very well with the determination by Kallivayalil et al. (2018).

A.20. Hydrus I

In the case of Hydrus I a spectroscopic catalog with continuous membership probabilities is published in the electronic table of Koposov et al. (2018). Of the 33 candidates, 30 remain after our cuts. One remaining member has a relatively low membership probability of 56%; if excluded, the motion changes by 0.1 and 0.4σ in RA and Dec, respectively. Kallivayalil et al. (2018) added some more photometric candidate members. Our motion agrees with their motion and also with the Simon (2018) values within about 1σ.

A.21. Leo I

For Leo I we use the catalog of Mateo et al. (2008), which does not contain membership probabilities or gravities. Therefore, we assign a probability on the basis of the l.o.s. heliocentric velocity. After our cuts we have 241 members, and thus the systemic motion is robust. We find good agreement with the Gaia Collaboration (2018b) determination.

A.22. Leo II

For this galaxy we use primarily the catalog of Spencer et al. (2017), which has binary probabilities. In addition, we complement the analysis with the catalog of Koch et al. (2009); to these stars we assign membership probabilities on the basis of l.o.s. heliocentric velocities. Combining both catalogs, 131 members pass all our cuts. Our motion agrees well with that of Gaia Collaboration (2018b).

A.23. Leo IV

For this galaxy we use the catalog of Simon & Geha (2007), which provides binary probabilities. Due to the galaxy distance and faintness, only three members remain after our cuts. A few faint stars without spectroscopy have a consistent motion to our measurement. Also, the background is low at the large distance of this galaxy. Due to the faintness of the three stars the error of the motion is still rather large and does not constrain the properties of the system well.

A.24. Leo V

We use the catalog of Walker et al. (2009a) with continuous membership probabilities. After our cuts we have five member stars, all with a high spectroscopic membership probability, which should imply they are robustly classified as members. In addition to the stars with spectroscopy, there is one faint star without spectroscopy with a consistent proper motion. Due to the distance and faintness of the system, the error on the systemic PM is rather large and does not constrain the orbital properties of the galaxy very well.

A.25. Phoenix I

We use the catalog of Kacharov et al. (2017), which provides binary probabilities. Because the galaxy is rather bright we have a relatively large number of members (71) even if the system is very distant. There are also more stars without spectroscopy that have proper motion consistent with our determination. The spectroscopy-based motion still does not constrain the properties of the galaxy very well.

A.26. Pisces II

We use the catalog of Kirby et al. (2015), which provides binary probabilities. Only two stars pass our cuts (9004 and 12924 in the nomenclature of Table 2 of Kirby et al. 2015) and there are no additional stars without spectroscopy that we can use to back up our determination. This is partly expected due to the faintness of the galaxy. Clearly, the systemic PM is less reliable than most of our determination, but its error is rather large.

A.27. Reticulum II

We use the catalog of Simon et al. (2015), giving binary membership classifications. After our cuts we have 27 member stars. Our motions agree very well with Simon (2018), Kallivayalil et al. (2018) and Pace & Li (2018), which agrees within 1.7 and 0.3σ, respectively with the motion of Massari & Helmi (2018).

A.28. Sagittarius I

For the Sagittarius dSph we use APOGEE DR14 data (Majewski et al. 2017), selecting all the stars within 30’ of the nominal galaxy center. To assign membership probabilities we use a selection based on heliocentric l.o.s. velocity and stellar gravity, but retaining only very bright giants; due to the location of Sagittarius close to the bulge and due to the selection function of APOGEE, nonmember giants were also observed. Sag I is one of the galaxies were the systemic PM changes more noticeably (by −1 and −0.6σ in RA and Dec) if we were to adopt a stricter membership cut; nonetheless, this is just a change of 0.01 mas yr⁻¹ or less. Our systemic PM differs from the determination in Gaia Collaboration (2018b) by 5σ (0.044 mas yr⁻¹ in μ_α). However, the difference is not much larger than the systematic error, which could play a role for this large galaxy. It is likely that the fact that Sag I is rather close also causes measurable proper motion gradients throughout the object, which could be present in varying amounts in samples with different spatial distributions. In order to achieve a higher accuracy, it might be necessary to consider a model accounting for possible such gradients.

A.29. Sculptor (I)

As for Fornax, we use the catalogs of Walker et al. (2009b) and Battaglia & Starkenburg (2012) and the same approach to member selection. We find good agreement with the motion of Gaia Collaboration (2018b).

A.30. Segue 1

We use the catalog of Simon et al. (2011) and the membership probabilities based on a Bayesian approach. Of the selected 13 stars, the star with the lowest probability has a p= 0.885, thus the adopted cut does not matter. The same is true if we use one of the other two sets of probabilities provided in the article. Our motion agrees within about 1σ with Simon (2018). It is different by 2.2 and 0.2σ (in RA and Dec) from the proper motion measurement of Fritz et al. (2018b), based on LBT data.

A.31. Segue 2

We use only the stars classified as certain members in Kirby et al. (2013). After our cuts, we are left with ten stars, which however do not form a clear clump in proper motion space. Our systemic PM agrees in RA with the motion of Simon (2018) but differs by 2σ in Dec, while it differs by 3.1 and 3σ (in RA and Dec) from the determination of Massari & Helmi (2018). In contrast to our work, the latter work also includes stars without spectroscopy. While adding these additional stars improves precision, it is not clear whether this is also true when uncertain membership is taken into account. We note that there is a stream at a similar distance to that of Segue 2 (Belokurov et al. 2009; Kirby et al. 2013), and that the use of projected distance from the galaxy center as a criterion to start isolating targets by Massari & Helmi (2018) is based on a simple cut. Because Segue 2 is at a close distance, pollution by other halo components is more important for this object than for most of the other systems in our sample. For an improved determination of the systemic PM, it might be necessary to wait for future Gaia data releases in which the proper motion errors are smaller and thus the member stars would form a tighter clump.

A.32. Sextans (I)

We use the catalog of Cicuéndez et al. (2018) with continuous probabilities. With our standard cuts we select 325 members. The proper motion changes by less than 0.1σ when we change the spectroscopic probability cut to P >  0.7. Our motion differs by 1.6 and 0.7σ (in RA and Dec) from the motion of Gaia Collaboration (2018b), which might happen by chance or due to additional systematic differences due to the large angular size of this galaxy.

A.33. Triangulum II

We use the catalog of Kirby et al. (2017), which provides binary memberships. After our cuts we are left with five members; the star classified as an uncertain member in the source paper is eliminated by our cuts. These stars do not form a clear clump in proper motion space, but due to the faintness of the stars and to the possible larger population due to Tri II low Galactic latitude, a large scatter is not surprising. Our motion agrees very well with the motion of Simon (2018).

A.34. Tucana II

We use the catalog of Walker et al. (2016). It provides nonbinary memberships. With our standard cuts we obtain 19 members stars. Most of those form a clump in proper motion space. Some of the others are more scattered, but since they are fainter they do not matter much for the overall motion. Our motion agrees within 1σ with the determination of Simon (2018), but differs by 0.9 and 3.0σ in RA and Dec with the determination by Kallivayalil et al. (2018). Our values agree nearly perfectly with those of Pace & Li (2018). This is evidence that the method of Kallivayalil et al. (2018) of adding photometric members is less reliable than the photometric method of Pace & Li (2018).

A.35. Tucana III

This satellite has a short tidal stream (Drlica-Wagner et al. 2015) that was covered spectroscopically by Li et al. (2018a). Along streams of globular clusters (we note that Tuc III has at most a dispersion of 2 km s⁻¹ (Li et al. 2018a; Simon et al. 2017), thus is dynamically very similar to a globular cluster), gradients in proper motions are smooth and show only a small jump around the main body; for Pal 5, which is closer in distance than Tuc III, this is less than 0.1 mas yr⁻¹. Such a contribution would be of the same order as systematic uncertainties for Gaia (Lindegren et al. 2018) over small areas; therefore, we decided to include stream stars. To avoid that a somewhat uneven coverage of the two sides of the stream could have some influence on the results, we fit the proper motion as a function of position and use here only the proper motion at the intercept, i.e., at the position of Tuc III. The systemic proper motion obtained in this way is consistent with the values in Simon (2018), who do not use stream stars, and in Pace & Li (2018), who use fewer stream stars than we do. Our results is also consistent with Kallivayalil et al. (2018), who obtain the proper motion in a more simple way.

A.36. Ursa Major I

We use the catalog of Martin et al. (2007) and, for stars not contained in that work, the catalog of Simon & Geha (2007). Both provide binary classifications. After our cuts we have 23 members. Our motion agrees very well with the determination of Simon (2018).

A.37. Ursa Major II

We use the catalogs of Martin et al. (2007) and, for stars not contained in that work, the catalog of Simon & Geha (2007). Both provide binary classifications. After our cuts we have 15 members. These stars form a clear clump in proper motion space supported by some stars without spectroscopy. Also in this case our motion agrees very well with the determination of Simon (2018).

A.38. Ursa Minor (I)

We use the catalog of Kirby et al. (2010), which supposedly only contains members. After our cuts we have 137 stars. The systemic PM agrees very well with the determination by Gaia Collaboration (2018b).

A.39. Willman 1

We use the catalog of Martin et al. (2007), which provides binary memberships. After our cuts we have seven member stars. The brightest of these stars form a group in proper motion space, which is supported by one star without spectroscopy. Our motion agrees within 1σ with the determination of Simon (2018), but our error is relatively large.

Appendix B: Monte Carlo simulations for errors

B.1. Forward MC simulations

We ran 2000 Monte Carlo simulations drawing μ_α^*, μ_δ, and the heliocentric l.o.s. velocities randomly from Gaussian distributions centered on the measured values and dispersions given by the respective errors. When extracting random values for μ_δ we consider the correlation between μ_α^* and μ_δ listed in Table 2. We also use the listed value for the systematic error. The randomly simulated kinematic and positional properties are then transformed into velocities in a Cartesian heliocentric (and then Galactocentric) reference system. From these values we then calculate orbital poles. As we tested in Monte Carlo simulations, the derived quantities mentioned above are not biased since they can take both positive and negative values.

B.2. Backward MC simulations

As shown in Fig. B.1, the process of normal forward MC simulations described above obtains biased errors and values for positive defined properties like ν_tan and ν_3D when errors are close to being as large as or larger than the observed value; this occurs because the error already affects the observed value and the standard MC simulations would add the error a second time. Therefore, we adopt a different approach for estimating ν_tan and ν_3D, and the corresponding errors. Essentially, we want to determine what combinations of proper motions and l.o.s. velocities would result in the observed values of ν_tan and ν_3D when they are convolved with their observational errors (see also van der Marel & Guhathakurta 2008 for an independent development of a similar methodology). To this end, we proceed in the following way. We draw values of μ_tot (i.e., $\sqrt{{\mu }_{\alpha ,\ast }^2+{\mu }_{\delta }^2}$) from a large grid of values, assuming a uniform distribution in μ_tot and in its direction. In addition, we draw most of the other properties as previously (e.g., l.o.s. velocity), with Gaussian errors. We then add the observed errors to μ_α, * and μ_δ and convert the proper motions and l.o.s. heliocentric velocities into ν_tan and ν_rad. At this point, we require that both velocities agree with their corresponding observed values within 0.05–0.8σ. We choose a tolerance of these sizes because for larger tolerances the retrieved values are bigger because numbers are somewhat more likely drawn on the larger side of the measured value. Our tolerance is larger for cases that have a more significant motion because biases are smaller in this case. We note that here we neglect the error in distance modulus to prevent the velocity selection from also selecting in distance. For the errors we use those obtained in the forward MC simulations. For all the galaxies we obtain at least 2000 accepted Monte Carlo values; this gives us a distribution of 2000 error-free proper motions that, when convolved with their errors (and associated with the l.o.s. velocity and distance of the galaxy), would produce unbiased tangential and 3D velocities. The distribution of 2000 simulated ν_rad, ν_tan, and ν_3D is then simply derived by transforming the 2000 error-free proper motions and Gaussian distribution l.o.s. velocities into Galactocentric velocities; in this step, we allow for distance errors. We use quantiles to give 1σ errors in ν_rad, ν_tan, and ν_3D, and use as a most-likely value the median of the distribution. For ν_rad the errors are nearly symmetric, thus we give a single error; for the other two we give asymmetric errors (see Table 2).

Fig. B.1. Comparison of errors and values for positive-defined quantities such as ν3D and νvtan, using forward and backward Monte Carlo simulations. The x-axis gives the ratio of observed velocity to the error given by the forward MC simulations, while the y-axis shows the ratio of the median velocity obtained from the MC simulations (forward or backward, see legend) to the observed value.

We compare in Fig. B.1 the difference between forward and backward simulations. The outcome mainly depends on the error in proper motion, which is closely related to the error in tangential velocity; of smaller importance are the errors in distance and l.o.s. velocity. It is clear that for ν_tan of less than about 4σ significance the tangential velocity is overestimated in the forward Monte Carlo simulations. The situation is similar for ν_3D, but since the l.o.s. velocity matters much more in this quantity, the relationship is less linear.

The bias in ν_tan and ν_3D also affects the estimate of orbital properties, such as apocenter and pericenters, because galaxies with a 3D velocity that is artificially too large would appear to be less bound than they really are. Therefore, we use the backward Monte Carlo simulation results (namely μ_α*, μ_δ, and V_LOS) also in Galpy, selecting for each satellite a random subsample of 500 cases. We then retrieve the median and the 1σ range in the distribution of peri-/apocenter and eccentricity from the 500 orbital integrations.

We note that because of the way the backward Monte Carlo simulations had to be set up, the orbit has a random orientation with respect to the disk, which can impact orbits with a small pericenter. Tuc III is the only satellite which is certainty affected; however, since its pericenter is just 3 kpc, other approximations, such as the use of a spherical bulge in Galpy, would also matter, making the used gravitational potential not fully realistic. All the other satellites most likely have pericenters of at least 14 kpc, such that the disk is less important than the positive bias of the forward Monte Carlo calculation.

All Tables

All Figures

	Fig. 1. Total velocity of all galaxies in the sample (points with error bars). The curves show the escape velocity for the two potentials used in this paper, i.e., MWPotential14 with a NFW halo of virial mass 0.8 × 1012 M⊙ (black line) and a more massive variant with virial mass 1.6 × 1012 M⊙ (red line).
In the text

	Fig. 2. Pericenter vs. apocenter properties of the galaxies in the sample with 3D velocity errors < 100 km s−1, as indicated by the labels. The top panel plot shows the results for a MW DM halo of virial mass equal to 1.6 × 1012 M⊙, and the lower panel for 0.8 × 1012 M⊙. We also show typical virial radii in both plots. Gru I is omitted in the lower plot because of the very large apocenter.
In the text

Fig. 3. All-sky view of orbital poles for the objects in the sample; the circles indicate the median of the 2000 Monte Carlo simulations, while the small points around each object plot the orbital poles from the individual simulations. The magenta circles contain 10% of the sky around the assumed VPOS pole, which is given as an “X” for the co-orbiting direction (orbital sense as most classical satellites, including the LMC and SMC), and a “+” for the opposite normal direction (counter-orbiting). Top panels: objects with Galactocentric distances between 0 and 50 kpc (left) and 50 and 100 kpc (right); bottom panels: objects with Galactocentric distances between 100 and 200 kpc (left) and in the putative Crater-Leo group (right).

In the text

	Fig. 4. Properties of Crater II. Top: orbital integration for Cra II in the MWpotential14 for the standard light halo (blue line) and the heavy DM halo (red line) plotted as the distance from the MW center as a function of time. Bottom: hill radius of a galaxy of the mass of Cra II (Caldwell et al. 2017) in two different MW potentials compared to the size and apocenter of Cra II.
In the text

	Fig. 5. Histograms of the ratio f = (dGC − rperi)/(rapo − rperi) for high-mass and low-mass MW models (upper and lower panels, respectively). The galaxy on the left of each histogram is near pericenter, while the galaxy on the right of each histogram is near apocenter.
In the text

Fig. 6. Alignment with the VPOS, measured as the ratio between the measured (θobs) and the predicted (θpred) angular offset of an orbital pole from the VPOS, plotted against the uncertainty (Δpole) in the direction of the orbital pole. Blue upward triangles are co-orbiting relative to the majority of classical satellites including the LMC; red downward triangles are counter-orbiting. Symbol size increases with Galactocentric distance.

In the text

	Fig. B.1. Comparison of errors and values for positive-defined quantities such as ν3D and νvtan, using forward and backward Monte Carlo simulations. The x-axis gives the ratio of observed velocity to the error given by the forward MC simulations, while the y-axis shows the ratio of the median velocity obtained from the MC simulations (forward or backward, see legend) to the observed value.
In the text